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CHAPTER  1 


INTRODUCTION 

1 , 1  Autonomous  Satellite  Navigation  History 

Self-contained  or  autonomous  navigation  of  spacecraft  was 
a  desired  capability  almost  at  the  very  beginning  of  space  flight. 
The  earliest  references  to  space  navigation  discuss  its  necessity 
for  manned  missions  and  interplanetary  travel  [Henry,  1963;  Gersten 
and  Schwarzbein,  1963],  but  unmanned  earth  satellites  have  continued 
to  be  tracked  and  controlled  by  worldwide  tracking  networks.  These 
networks,  set  up  by  both  the  United  States  Air  Force  and  the 
National  Aeronautics  and  Space  Administration  (naSA)  are  complex, 
expensive,  redundant  and  require  large  operation  and  maintenance 
budgets . 

Studies  of  artificial  satellite  autonomous  navigation 
sensors  and  techniques  have  been  performed  since  the  late  I960 
decade  by  Brogan  and  Lemay  [ 1 968 J ,  Gura,  et  al .  [ 1 97 1  ]  and  Lemay,  et 
al.  [1973],  but  sensor  development  has  lagged  far  behind  the 
analyses  and  only  recently  have  serious  moves  beer,  taken  to  build 
and  test  sensors  that  will  enable  transfer  of  the  navigation 
function  from  the  earth-based  system  to  each  active  soacecraft. 
NASA,  driven  primarily  by  cost  considerations,  is  planning  to  jse 
the  Tracking  and  Data  Relay  Satellite  System  (TDRSS)  as  an  orbital 
tracking  station,  with  data  processing  still  to  be  performed  on  the 
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ground.  The  Department  of  Defense  (DOD),  however,  is  more  concerned 
with  vulnerability,  as  stated  by  Robert  S.  Cooper,  head  of  the 
Defense  Advanced  Research  Projects  Agency  [Aviation  Week.  1982],  and 
with  the  overseas  ground  stations  being  the  most  vulnerable  links  in 
the  tracking  and  control  system,  DOD  is  funding  spaceborne  sensor 
systems  of  varying  degrees  of  autonomy  to  enable  the  onboard 
performance  of  the  navigation  function. 

Previous  autonomous  satellite  navigation  investigations 
involved  several  different  types  of  sensors.  Lemay,  et  al.  [ 1 973 ] 
thoroughly  investigated  the  use  of  both  known  and  unknown  landmark 
trackers,  horizon  scanners,  satellite-to-satellite  measurements  of 
angles,  range  and  range-rate,  star-horizon  sensors  and  space  sextant 
measurements.  Their  investigation,  which  was  based  upon  state-of- 
the-art  sensor  precision  in  1973,  indicated  that  landmark  trackers 
had  the  potential  for  yielding  the  best  navigation  performance. 
Development  of  landmark  trackers,  however,  never  achieved  the 
potential  expected  of  them,  and  optical  star  trackers,  along  with 
the  Global  Positioning  System  (GPS),  are  the  onboard  navigation 
systems  currently  under  development. 

It  is  interesting  to  note  that  USSR  interest  in  autonomous 
navigation  started  somewhat  parallel  to  but  behind  the  United 
States.  Tne  paper  by  Zybin  [1969]  proposes  using  star-planet 
observations  in  a  deterministic  orbit  determination  scheme  similar 
to  Gersten  [1963].  On-orbit  testing  of  space  sextants  occurred 
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early  in  both  manned  space  programs.  Experiments  were  performed  on 
both  the  Gemini  [Ballentine,  196?]  and  Skylab  [Walsh  and  Ferrguson, 
1975]  spacecraft,  as  well  as  in  the  Soviet  manned  spacecraft 
[Nikoloev,  et  al.,  1975]. 

1 .2  Air  Force  Spacecraft  Navigation  Requirements 

The  new  sensors  under  development  promise  navigation 
accuracies  that  may  be  competitive  with  ground-based  systems  in 
meeting  most  current  and  projected  Air  Force  spacecraft  position 
;nowledge  requirements.  Discussions  with  Air  Force  System  Program 
Offices  (sPO's)  during  the  summer  of  1981  led  to  a  list  of  accuracy 
requirements  whose  one-sigma  values  spanned  the  range  from  less  than 
10  m  to  ore  than  37  km  [Tapley  and  Ferguson,  1983], 

To  meet  these  requirements,  some  of  the  program  otfices 
are  investigating  the  use  of  current  sensor  data,  while  others  will 
require  new  sensors  of  the  type  being  developed.  The  Defense 
Meteorological  Support  Program  (dmsp)  is  tyoical  of  those 
investigating  their  current  sensor  capabilities.  The  DMSP  Primary 
Attitude  Determination  System  (PADS)  consists  of  a  fixed  star 
tracker  pair,  a  sun  sensor  and  earth  horizon  sensors.  When  all 
sensors  are  In  operation,  the  attitude  error  sigma  is  36  arcsec  and 
the  system  outputs  attitude  errors  and  a  unit  vector  pointing  to  the 
center  of  the  earth.  This  unit  vector  can  be  utilized  by  a 
navigation  filter  to  perform  the  autonomous  navigation  function. 
Current  sensor  data  studies  have  a  goal  of  1  nm  accuracy  with 
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atmospheric  radiance,  bia3  modeling  and  filter  design  being  the 
areas  needing  improvements.  Since  the  program  navigation 
requirement  is  .5  nm  (930  meters),  improvements  to  the  sensors  or 
new  horizon  sensors  are  necessary.  Other  programs,  whose  accuracy 
requirements  are  in  the  100  m-1  km  range,  are  also  candidates  for 
the  optical  star  trackers  being  built  and  tested  under  Space 
Division  auspices. 

Some  satellite  programs,  however,  require  position 
knowledge  to  50  m  or  less  and  will  thus  necessitate  the  development 
of  extremely  precise  onboard  measurement  equipment.  The  GPS  program 
is  not  only  one  of  those  that  requires  very  accurate  navigation 
information  but  is  the  only  navigation  system  proposed  whose 
position  accuracies  meet  the  requirements  of  other  high-precision 
users.  It's  main  limitation  from  an  autonomous  system  viewpoint  is 
that  it  is  dependent  upon  a  global  tracking  system.  User  satellites 
are  thus  dependent  upon  a  navigation  system  that  is  still  vulnerable 
to  ground  system  failures  and  outages. 

The  GPS  Joint  Program  Office  (JPO)  has  performed  limited 
studies  on  two  autonomous  navigation  schemes,  one  using  range  and 
integrated  doppler  measurements  from  other  GPS  vehicles  via  the 
cross-link  antenna  and  the  other  involving  precise  horizon  data  from 
a  new  optical  sensor.  The  GPS  satellite  relative  range  study  by  Liu 
[ 1 981 j  indicates  that  although  in-track  errors  grow  to  30  meters, 
10  meter  accuracy  can  be  maintained  in  the  user  line-of-sight 
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direction  for  Phase  I  vehicles  involved  in  active  satellite-to- 
satellite  cracking  for  I1)  days  following  one  day  of  ground  tracking. 
Earth  geopotential  resonance  terms  are  the  largest  error  sources, 
and  improved  values  for  C^,  sw,  C32  and  S32  are  required.  The 
study  did  not  address  the  Phase  III  constellation,  long-term 
stability,  clock  variation  or  stability  of  an  onboard  solution 
process. 

The  optical  sensor  proposed  for  GPS  measures  the  earth's 
horizon  by  atmospheric  dispersion  of  star  images,  and  is  named  SHAD 
for  Stellar  Horizon  by  Atmospheric  Dispersion.  The  current  sensor 
design  is  intended  for  mid-course  missile  guidance  and  is  expected 
to  result  in  a  position  error  sigma  of  65  m  [Quell,  1981  ],  so  it 
would  require  approximately  an  order  of  magnitude  improvement  to  be 
used  by  GPS.  Earth  atmospheric  density  modeling,  as  discussed  in 
Section  2.8,  is  a  major  factor  limiting  the  accuracy  achievable  by 
SHAD. 

1 -3  Purpose  of  the  Study 

While  some  Air  Force  satellite  SPO's  are  actively 
investigating  the  capability  of  current  and  future  onboard  sensors, 
several  others  are  either  independently  soliciting  proposals  for  new 
sensors  or  are  waiting  for  sensors  to  be  developed  to  the 
operational  stage  before  making  a  decision  on  the  route  they  will 
take  to  satisfy  autonomous  navigation  requirements.  There  is  a  need 
to  assess  the  capabilities  of  sensors  now  under  development,  match 
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them  to  these  program  requirements  and  identify  needed  improvements 
so  that  decisions  can  be  made  as  to  the  direction  future  sensor 
development  should  take.  The  purpose  of  this  study,  then,  is  to 
assess  navigation  accuracies  using  these  new  sensors  so  they  may  be 
matched  to  existing  and  planned  satellite  missions.  Since  some 
users  have  requirements  that  can  only  be  met  by  GPS,  an 
investigation  of  improved  GPS  autonomy  is  also  undertaken. 

The  report  consists  of  a  description  of  the  new 
technology,  a  comparison  of  the  two  main  types  of  optical  sensors 
under  development:  the  Space  Sextant  and  the  matrix  star  sensor, 
and  an  analysis  of  GPS  autonomous  navigation  using  satellite-to- 
satellite  range  and  integrated  doppler  measurements.  The  GPS 
section  contains  a  relative  geometry  description,  a  description  of 
relevant  error  sources,  ephemeris  model  selection  using  consider 
analysis,  orbit  and  clock  simulation  descriptions  and  analysis  of  a 
proposed  local  estimation  algorithm.  Conclusions  concerning  the 
various  sensors  are  drawn  and  recommendations  are  made. 

1 . M  Description  of  Sensors 

Five  of  the  sensors  now  under  development  yield  data  which 
hold  high  promise  for  current  autonomous  navigation  applications. 
These  include  the  Space  Sextant  Autonomous  Navigation  and  Attitude 
Reference  System  (SS/ANARS),  the  Multimlssion  Attitude  Determination 
and  Autonomous  Navigation  System  (MADAN)  and  the  Digistar  and 
STELLAR  star  sensors.  The  measurement  from  each  of  these  sensors  is 
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based  on  the  sensor's  ability  to  measure  accurately  the  angle 
between  a  star  and  a  near  celestial  body  or  the  angle  between  two 
stars.  One  additional  3ensor,  the  GPSPAC,  uses  range  and/or 
integrated  doppler  measurements  from  the  Global  Positioning  System 
satellites  as  data  for  the  navigation  function.  The  advertised 
characteristics  as  determined  by  the  design  specifications  for  these 
sensors  are  given  in  Table  1.1.  Operational  characteristics  of  each 
sensor  are  described  in  further  detail  in  the  subsequent  discussion. 


TABLE  1.1.  ADVERTISED  SYSTEM  CHARACTERISTICS 


IZE £ 


1 c  Sensor 

Weight 

Power 

Precision 

Operational 

Date 


Space  Sextant 

MADAN 

Digistar 

STELLAR 

GPSPAC 


.5  arcsec 
2  arcsec 
.5-2  arcsec 
1-30  arcsec 
Pos:  18  m 
Att:  .02°-. 


65  lb3 
50  lbs 
30-60  lbs 
40  lbs 
43  lbs 
0  (71-2143 


50  watts 
50  watts 
30-60  watts 
40  watts 
45  watts 
arcsec) 


1985 

1987 

1985  (2  arcsec J 
1987 


1.4.1  Space  Sextant 

The  Space  Sextant  Autonomous  Navigation  and  Attitude 
Reference  System  is  being  developed  by  Martin-Marietta.  A  flight 
demonstration  model  has  been  built  and  ground-tested,  and  an 
operational  version  may  be  ready  by  1985.  The  test  model  weighs 
approximately  220  lbs  and  consumes  280  watts,  but  the  operational 
model  proposed  for  the  Mini-HALO  program  is  advertised  to  weigh  65 
lbs  and  consume  50  watts  of  power  [Martin-Marietta,  1 980 ] . 
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The  sext&it  is  composed  of  two  Cassegrain  tracking 
telescopes  mounted  on  a  3  degree- of-freedom  inertial  platform  (Fig. 
1.1 ).  In  che  navigation  mode,  one  telescope  tracks  the  bright  limb 
oi  the  moon,  while  the  other  tracks  stars  visible  to  the  system.  A 
timing  wheel  located  between  the  telescopes  rotates  a  prism  at  9  rps 
such  that  optical  signals  are  injected  into  each  telescope  parallel 
to  the  received  starlight.  The  angle  between  tne  two  lines  of  sight 
is  then  determined  to  <  .5  arcsec  by  measuring  the  time  interval 
between  the  optical  signal  reception  at  each  telescope's  detector. 

1 . .  2  MAPAN 

The  Multimission  Attitude  Determination  and  Autonomous 
Navigation  system  is  a  solid-state  matrix  star  sensor  being- 
developed  by  TRW  (Fig.  1.2).  The  heart  of  the  sensor  is  a  matrix 
charge-coupled  device  (CCD)  developed  by  Hughes.  The  matrix 
contains  four  arrays  of  324  x  324  elements  or  pixels,  each  1  mil  x  1 
mil  (25.4  um  x  25.4  pm)  [TRW,  1979].  A  Schmidt-Ca3segrain 
reflecting  telescope  with  a  7.1°  x  7.1°  field  of  view  produces  an 
intentionally  defocused  image  on  the  array  and  a  sensor  data 
processor  determines  the  centroid  of  the  image  to  approximately  5 
percent  of  a  pixel  width  with  respect  to  the  sensor  line  of  sight. 
Since  each  pixel  subtends  39.4  sec,  a  5  percent  error  gives  1o  <  2 
arcsec.  Two  such  sensors  can  be  used  to  determine  spacecraft 
attitude,  but  the  star  sensor,  along  with  its  data  processor,  is  not 
capable  of  autonomous  navigation  without  an  earth  horizon  sensor. 


It  may  be  possible  to  use  the  star  sensor  as  an  earth  horizon  sensor 
by  measuring  a  star  position  as  it  is  refracted  by  the  atmosphere, 
but  the  v.urrent  design  does  not  include  the  necessary  software.  The 
use  of  this  sensor  for  horizon  determination  is  discussed  further  in 
Section  2.4.1. 

A  MADAN  test  model  is  being  fabricated  by  TRW,  with  bench 
test  results  expected  in  1983.  An  operational  version  could  be 
flown  by  1987. 

1.4.3  Dlglstar 

Digistar  is  another  solid-state  matrix  star  sensor  being 
built  by  Ball  Aerospace  Systems  Division  (Fig.  1.3).  It  employs  a 
256  x  256  pixel  charge  injection  device  (CID)  developed  by  the 
General  Electric  Corporation.  Each  pixel  element  is  20  x  20  pm. 
The  refracting  telescope  produces  an  image  on  the  focal  plane 
matrix,  and  an  interpolation  scheme  is  used  to  arrive  at  an  image 
centroid  with  a  precision  of  less  than  1  percent  of  pixel  size  in  a 
field  of  view  of  2.93°  x  2.93°.  Since  each  pixel  subtends  41 
arcsec,  the  resulting  star  position  precision  is  expected  to  be 
about  .4  arcsec.  Testing  of  a  breadboard  model  currently 
demonstrates  a  1o  tracking  error  of  .8  arcsec  [Ball  Aerospace, 
198l].  According  to  a  Ball  representative,  an  operational  model 
could  fly  in  1985. 
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1.4  4  STELLAR 

The  Jer  Propulsion  Laboratory  (JPL)  has  been  developing 
CCD-based  star  sensors  since  1974  [Salomon,  1981  ].  The  project  is 
named  STELLAR  for  Star  Tracker  for  Economical  Long  Life  Attitude 
Reference.  To  date,  tv;o  large  array  trackers  have  been  developed. 
The  Video  Inertial  Pointing  [VIP)  tracker  employs  a  Fairchild- 
produced  CCD  having  a  190  x  244  element  array  with  a  1.9°  x  2.5° 
f ield-of-view,  36  arcsec  pixels  and  2.2  arcsec  resolution.  This 
instrument  flew  on  a  balloon  payload  in  June  1979. 

The  Extended  Life  Attitude  Control  System  (ELACS)  tracker 
uses  a  38O  x  488  element  CCD  built  by  Fairchild  and  has  a  )0°  x  32° 
f ield-of-view.  The  resolution  is  12  arcsec  along  the  short  axis  and 
30  arcsec  along  the  long  axis.  Even  though  it  is  less  accurate  than 
the  VIP  tracker,  it  allows  commandable  f ields-of-view,  so  several 
operating  modes,  such  as  star  field  mapping,  star  acquisition  and 
star  tracking,  are  available  from  one  instrument. 

The  JPL  goal  is  to  have  three  instruments:  a  wide-field 
star  tracker,  like  ELACS,  a  general  purpose  tracker  and  an 
instrument  pointing  sensor  with  accuracies  varying  from  30  to  1 
arcsec.  According  to  Salomon  [1981],  each  tracker  would  weigh 
between  8  and  13  lb  and  draw  10  watts  of  power. 

1.4.5  GPSPAC 

A  spaceborne  Global  Position  System  receiver  (GPSPAC)  is 
under  development  by  Magnavox,  among  others.  Utilizing  the  GPS 
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ranging  signal,  position  errors  of  18  m  (lo)  are  expected  for  user 
satellites  below  GPS  altitude  (20197  km)  in  Phase  III  operation. 
This  system  is  dependent  upon  the  GPS  spacecraft  launch  schedule  but 
is  expected  to  be  operational  in  1987. 

User  satellite  attitude  can  be  determined  from  GPS-derived 
interferometric  information  obtained  by  using  two  widely  separated 
antennas  and  a  suitably  modified  receiver.  This  information  is 
obtained  separately  from  the  navigation  data  and  requires  additional 
receiver  design.  Attitude  determination  accuracy  using  this 
technique  is  only  on  the  order  of  .02°  to  .6°  [Ellis  and  Creswell, 
1978]. 

1 .5  Sensor  Comparison 

The  relative  advantages  and  disadvantages  of  each  sensor, 
from  the  operational  point  of  view,  can  be  summarized  as  follows: 

A.  Space  Sextant 
Advantages: 

1 .  High-accuracy  angle  measurement 

2.  Self-gimbaled 

3.  Early  operational  date 
Disadvantages: 

1.  Mechanical  gimbals,  possibly  reducing  reliability 

2.  Moon-star  angle  is  less  sensitive  to  orbit  dynamics  than 
an  earth-scar  angle. 

B.  Solid-state  Matrix  Sensor  (MADA rt,  Digistar  and  STELLAR) 
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Advantaged. 

1 .  Sinai1  ,  30I i d  state 

2.  Modular 

i.  t.an  view  several  stars  simultaneously 
Disadvantages: 

1.  Fixed  to  spacecraft  with  no  automatic  scan  capability 

2.  No  existing  comparable  earth  horizon  sensor 

3.  Later  operational  date  for  a  navigation  system 
C.  GPS  Receiver 

Advantages:  High  positional  accuracy 
Disadvantages: 

1.  Dependent  on  ground  tracking  system 

2.  Low  attitude  precision 


CHAPTER  2 


OPTICAL  SENSOR  NAVIGATION  ACCURACIES 

2 . 1  Introduction 

Since  the  various  instruments  described  in  the  previous 
sections  will  provide  different  navigation  information,  conversions 
must  be  made  from  sensor  precision  to  spacecraft  position  error  to 
compare  the  navigation  accuracies  of  the  three  types  of  instruments. 
For  the  GPSPAC,  simulations  of  the  ranging  system  errors  coupled 
with  GPS  ephemeris  and  clock  errors  produce  a  user  position  error  of 
less  than  18  meters  during  periods  when  at  least  four  satellites 
with  acceptable  geometrical  displacement  are  visible  [Fuchs,  et  al . , 
1977].  Since  earth-orbiting  satellites  will  continuously  see  at 
least  four  GPS  vehicles  in  Phase  III,  this  accuracy  is  assumed  for 
all  users.  Semi-autonomous  maintenance  of  GPS  ephemerides  is 
analyzed  in  later  chapters. 

The  Space  Sextant  has  been  analyzed  in  detail  by  the 
Martin-Marietta  Corporation  [ 1 975  ] ,  and  a  navigation  accuracy  of 
300  meters  is  predicted  for  a  .5  arcsec  sensor.  The  matrix  sensors 
have  not  been  analyzed  in  detail  nor  have  the  comparisons  between 
them  and  the  Space  Sextant  been  made.  Such  a  comparison  will  be  a 
primary  objective  of  this  study.  Since  the  expected  navigation 
accuracy  of  the  Space  Sextant  is  known,  the  information  content  of 
the  star-moon  measurement  can  be  compared  with  that  of  the  star- 
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earth  horizon  observations  produced  by  the  matrix  sensors,  and  the 
resulting  relative  accuracies  will  be  indicative  of  matrix  sensor 
performance. 

2.2  Optical  Satellite  Navigation  Covariance  Program 

To  compare  the  navigation  accuracies  of  the  two  types  of 
optical  sensors  involved  in  this  study,  a  computer  program  for 
performing  covariance  propagation  and  analysis  was  written.  This 
program  assumes  only  zero-mean  Gaussian  measurement  errors  are 
present  and  propagate  the  satellite  state  error  covariance  using  the 
following  Extended  Kalman  Filter  (EKF)  update  scheme.  The  nonlinear 
model  for  the  system  is  represented  by 

X  -  F(x(t),t)  (2.0 

where  the  state  vector,  x(t),  has  components  measured  in  an  inertial 
frame: 

x,(t) 
x2(t) 
t) 
t) 
x5(t) 

X6(0 


x(t) 


x3( 

X„( 


F(t) 

v(t) 


(2.2) 
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and  where  the  force  function,  F(x(t),t],  is  defined  as: 


x4<t) 

I 


I 


1  r(0 

FvX(t).t)  =  J  .. 

[v(t)J 


x6(t) 

-px  (t) 

— V  + vc) 

r 

-Px?(t) 

- — +Py(t> 

r  J 

-Px  (t) 

- V"+  pz(t) 


(2.3) 


where 

o  2 

y  «  gravitation  parameter  of  earth  »  398603.2  km  /sec  , 

r  »  the  magnitude  of  the  satellite  position  vector, 

‘  (*?  *  -2  • 

P(t)  is  any  perturbing  force. 

In  this  analysis,  P(t)  «=  0  ,  and  the  resulting  orbit  is  two-body; 
however,  the  program  can  simulate  any  desired  forces  acting  on  the 
spacecraft  by  changing  the  derivative  subroutine. 

Between  planned  observation  times,  the  state  is  propagated 
by  either  closed-form  analytic  integration  for  the  two-body  case  or 
by  a  Runge-Kutta  (RK4)  or  Adams  predictor-corrector  numerical 
integrator.  The  state  error  covariance  matrix,  P,  is  given  at  the 
start  of  the  analysis  as  PQ  at  t  and  is  propagated  between 


observation  times  by 
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F(t  +  At  J  -  <t(  t  +  At,t )  P(t)  4>T(  t  +  At ,  t ) 


(2. A) 


where  <t> ( t  *■  At,t)  ,  t lie  state  transition  matrix,  is  approximated 
from  a  second-order  Taylor  series  solution  to 


4>(t,t)  -  A ( t )  «  (:,t);  4>(t,t)  -  I 


(2.5) 


where 


3X(t) 

The  approximate  solution  for  <!>(t  +  At,t)  is  given  by 


Murata  [ 1 982 ]  as 


$(t  +  At , t )  -  I  +  A(t ) At  +  [ A ( t )  +  A2 ( t )  ]  At' 


(2.6) 


where 


A ( t )  ■  ~  AUJ  (2.7) 

At 

At  each  simulated  observation  time,  tk  ,  the  covariance 
matrix  is  updated  by  the  contribution  of  the  type  of  observation 
assumed  to  be  available  at  that  time. 


2.2.1  Measurement  Model 

At  each  measurement  epoch,  the  covariance  matrix  is 
updated,  and  symmetry  is  enforced  by  the  following  algorithms: 


p  -  (1  -  kh)f 


(2.8) 


P  :*■  i  (p  +  PT) 


(2.9) 


Figure  2.1.  Spacecraft  Sensor  Orientation 
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are  related  to  the  sensor-orbit  geometry  and  are  usually  aligned 
with  the  largest  axis  along  the  T  unit  vector  of  the  RTN  system. 
The  velocity  error  ellipsoid  given  by 


defines  an  instantaneous  inertial  velocity  error  located  at  the 
satellite  but  does  not  include  terms  associated  with  the  angular 
velocity  and  acceleration  of  the  RTN  coordinate  system  itself. 

The  measurement  sensitivity  matrix,  H  ,  relates  the 
difference  between  the  measured  and  computed  observations, 
y  -  Y-Y*  ,  using  an  a  priori  estimate  of  the  trajectory,  to  the 
er-or  in  the  state  x  =  X  -  X*.  The  covariance  matrix  P  represents 
the  uncertainty  in  the  estimate  of  x.  The  structure  of  the 
measurement  se..„:,ivity  matrix,  H  ,  will  depend  on  the  individual 
observation  types.  For  the  cases  considered,  H  is  given  by: 
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3  measurement  (tk) 


ax. 


(2.14) 


2. 2. 1.1  Star-Horizon  Geometry 

For  a  star-horizon  observation,  Fig.  2.2  shows  the 
geometry,  with  vectors  and  angles  defined  as  follows.  Let 


(2.15) 


and 


rh  -  ray  tangent  aititude(RTA)  +  r 

§  =  unit  vector  to  star,  defined  by  line-of-sight  of  sensor 

in  inertial  coordinates 


Then, 


ob  «  b  -  a 


(2.16) 


(2.17) 

(2.18) 


and  the  first  three  components  of  the  measurement  matrix  are 
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while  the  second  three  components  are 


"*,5,6...  -»  (2'2°) 

2. 2. 1.2  Horizon  Sensor  Geometry 

For  a  horizon  observation,  the  simulated  star  direction  is 
assumed  to  be  in  the  plane  defined  by  the  center  of  the  earth,  the 
spacecraft  and  the  lxne-of-sight  of  the  same  star  sensor  simulated 
in  the  star-horizon  observation.  The  unit  vector  §  is  replaced  by 
a  unit  vector  pointing  to  the  intersection  of  this  plane  and  the 
surface  of  the  earth's  atmosphere  as  given  by  an  input  tangent 
height.  The  new  value  for  §  is 


sin(b) 


sin(a)  §  -  sin(ob)  — 
r 


(2.21) 


and  the  partial  derivative  matrix,  H  ,  is  the  same  as  aoove. 


2. 2. 1.3  Star-Moon  Geometry 

A  star-moon  observation  is  defined  as  shown  in  Fig.  2.3, 
with  the  observation  being  the  acute  angle  between  the  lunar  limb 


nearest  a  star  and  that  star. 


ob 


cos 


-1 


§  •  P„ 


-  sin 


and 


(2,22) 


H1,2,3 


1 


(2.23) 


H4,5,6  “  0 


(2.24) 


2 • 3  Star-Moon  Sensor  Performance 

The  space  sextant  uses  star-moon  angles  for  input  co  a 
navigation  filter  because  the  earth  horizon  cannot  be  accurately 
determined  by  current  infra-red  and  visible  sensor  technologies. 
Since  the  sextant  telescopes  can  gimbal  freely  with  respect  to  each 
other  and  with  respect  to  the  spacecraft,  one  can  normally  track  the 
bright  lunar  limb  while  the  other  locks  on  stars  visible  to  the 
system.  This  allows  a  more  precise  angle  measurement  than  is 
possible  with  a  current  earth  horizon  sensor.  In  the  situations  of 
lunar  occultation  by  the  earth  or  sun-moon  interference,  the  sensor 
tracks  the  earth  horizon  as  a  temporary  replacement  for  the  moon. 
The  sensor's  independence  of  motion  gives  it  the  advantage  of  being 
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abie  to  track  known  targets,  even  when  the  spacecraft  is  unstable  in 
attituae  and  also  frees  the  navigation  solution  from  attitude 
errors . 

The  space  sextant  type  of  observation  was  simulated  by 
computing  the  angla  between  unit  vectors  to  the  nearest  lunar  limb 
and  stars  at  ±45°  from  the  spacecraft  T  unit  vector  in  the  local 
horizon  plane.  As  shown  in  Figure  2.4,  the  resulting  navigation 
error  exhibits  a  large  twice-per-orbit  periodic  effect.  This  is  due 
to  the  orbit-moon  geometry  and,  for  various  initial  moon-ascending 
node  angles,  exhibits  varying  amplitudes.  The  case  with  minimum 
periodic  amplitude  (ft  =  73.4°)  was  then  U3ed  for  the  star-moon 
accuracy  curves  in  the  main  report  to  maintain  as  much  consistency 
as  possible  in  the  error  vs.  period  plots. 

This  twice-per-orbit  fluctuation  is  due  to  the  particular 
simulation  geometry  used  —  two  fixed  star  directions  relative  to 
the  spacecraft.  To  more  closely  simulate  the  Space  Sextant  scanning 
scheme,  as  described  by  Martin-Marietta  [ 1 98 1  ] ,  another  mode  was 
programmed  in  which  the  star  unit  vector  was  changed  randomly 
between  0°  and  360°  in  right  ascension  and  between  -90°  and  +90°  in 
declination  for  each  consecutive  simulated  measurement.  The  results 
(Fig.  2.5)  show  a  random  pattern  in  position  error  after  convergence 
of  the  filter.  This  random  behavior  would  adversely  affect  the 
consistency  of  results  when  various  period  orbits  are  considered,  so 
the  original  sensor  configuration  with  ft  *=  73.4°  was  retained  in  the 


analysis. 


2. 4  Star-Horizon  Sensor  Performance 

The  solid-state  modular  sensors  (MADAN,  Digistar  and 
STELLAR)  are  designed  to  be  rigidly  fixed  to  the  spacecraft  and, 
therefore,  cannot  observe  the  lunar  limb  except  when  it  happens  to 
fall  in  the  field  of  view.  Optical  celestial  navigation  from  a 
spacecraft,  however,  requires  the  measurement  of  angles  between  a 
star  and  a  body  that  appears  to  move  as  the  vehicle  moves  in  its 
orbit.  This  rules  out  star-star  angles  except  for  attitude 
information.  One  sensor  must  observe  a  near  celestial  body  and,  for 
earth-oriented  spacecraft,  the  earth  horizon  is  the  only  feasible 
target  for  fixed  sensors.  Also,  since  the  earth  horizon  appears  to 
move  360°  during  each  orbit  while  the  moon  appears  to  move 
2  tan  1(r/R[nj  (12.6  deg  for  sync  orbit),  the  earth-star  angle  is 
much  more  sensitive  to  orbit  dynamics.  Furthermore,  since  the  earth 
horizon  is  normally  closer  to  the  spacecraft  than  the  moon  horizon, 
the  earth-horizon  based  measurement  should  give  a  higher  navigation 
accuracy  for  angle  measurements  of  comparable  precision. 

For  near-term  conventional  horizon  sensors,  the  state  of 
the  art  is  about  .02°  [Fowler,  1981  ].  Figure  2.6  shows  the  results 
obtained  in  the  covariance  analysis  for  both  star-earth  and  star- 
moon  sensors.  The  covariance  analysis  program  propagates  the  state 
covariance  matrix  and  updates  this  matrix  with  information  obtained 
from  each  type  of  observation  at  fixed-time  intervals.  The 


resulting  FtSS  position  sigma  from  the  diagonal  terms  of  the 
covariance  is  plotted  versus  time  expressed  in  orbits,  with  the 
solid  line  representing  a  one-sigma  error  for  star-horizon 
measurements  obtained  from  a  2  arcsec  star  sensor  and  a  .02  deg 
horizon  sensor.  The  dashed  line  represents  the  one-sigma  error  for 
a  Space  Sextant  measuring  star-moon  limb  angles  to  a  .5  arcsec 
precision.  For  the  star-earth  measurement,  the  earth's  horizon  is 
sensed,  and  for  a  star-moon  angle,  the  moon’s  bright  limb  is  sensed. 
Two  horizon  sensors  with  .02°  accuracy  capability,  when  used  with 
two  fixed  star  sensors  with  a  precision  of  2  arcsec,  would  provide  a 
navigation  accuracy  of  2500  meters  for  a  satellite  in  a  12-hour 
circular  orbit.  As  seen  in  Fig.  2.6,  this  is  more  than  a  factor  of 
10  worse  than  that  achievable  by  a  space  sextant  with  a  precision  of 
.5  arcsec.  Thus,  to  take  advantage  of  the  improved  measurement- 
orbit  geometry,  the  solid  state  star  sensors  must  be  coupled  with 
horizon  sensors  of  at  least  an  order  of  magnitude  improvement: 
0.002°  or  7  arcsec. 

2.11.1  Star  Refraction  Measurement  of  Earth  Horizon 

Fortunately,  the  new  star  sensors  themselves  offer  a  means 
by  which  earth-horizon  sensing  may  be  improved.  As  starlight  passes 
through  the  atmosphere,  it  is  refracted  and  dispersed.  The  angles 
of  refraction  and  dispersion  depend  upon  ray  tangent  altitude  (RTA), 
the  point  where  a  starlight  ray  is  nearest  to  the  earth's  surface  as 
shown  in  Figure  2.7,  and  atmospheric  conditions,  but  the  measurement 
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Figure  2.6.  Navigation  Error  for  12h  Circular  Orbit;  Earth- 
Horizon  Sensor  Precision  «  0.02° 
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of  either  of  these  angles  offer  the  possibility  of  a  more  accurately 
defined  horizon.  The  Office  of  Naval  Research  is  presently 
investigating  the  feasibility  of  using  the  dispersion  sensor  SHAD, 
for  mid-course  ICBM  guidance,  with  an  expected  1o  error  of  65  meters 
[Quell,  1981  ],  while  Chambers  [ 1 981  ]  of  the  Aerospace  Corporation 
proposes  measuring  refraction  using  the  unique  capabilities  of 
MADAN.  Since  this  sensor  has  a  wide  field  of  view  (7.1°  square),  a 
minimum  of  two  stars  can  bfe  observed  continuously.'  By  measuring  the 
angle  between  two  stars  before  and  during  the  time  one  of  them  is 
refracted  by  the  atmosphere  due  to  spacecraft  motion,  the  refraction 
angle  inferred  from  this  changing  geometry  defines  the  height  of 
starlight  tangency  to  the  atmosphere  (RTA).  According  to  Chambers, 
a  3  arcsec  knowledge  of  this  refraction  angle,  obtained  by  looking 
at  two  stars  to  a  precision  of  2  arcsec  each,  gives  a  tangent  height 
uncertainty  of  150  meters  for  a  tangent  height  of  25  km.  Since  this 
uncertainty  depends  upon  a  priori  knowledge  of  the  properties  of  the 
atmosphere  and  sensor  and  not  on  orbit  altitude,  the  resulting 
apparent  horizon  determination  precision  varies  from  16.7  arcsec  for 
a  90-minute  orbit  to  .72  arcsec  for  a  geosynchronous  orbit. 
Figure  2.8  from  Chambers  [ 1 981 ]  shows  that  the  refraction  curve  has 
a  fairly  consistent  slope  for  all  models  considered  such  that  a 
given  refraction  error  produces  the  same  tangent  height  error  for 
each  curve.  In  addition,  this  tangent  height  error  is  only  a 
function  of  refraction  error,  and  refraction  error  is  a  function  of 
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the  optical  measurement  device,  not  orbit  altitude. 

When  the  solid-state  sensors  are  used  as  horizon  sensors, 
they  can  either  be  combined  with  other  similar  devices  to  provide 
star-horizon  data  or  they  can  be  employed  alone,  with  a  priori 
knowledge  of  the  angular  position  of  the  stars  being  observed,  to 
give  the  angle  between  the  horizon  and  the  center  of  the  earth. 
Fig.  2.9  shows  the  navigation  accuracies  expected  from  two  star- 
horizon  sensors,  two  horizon  sensors  and  two  star-moon  sensors  in  a 
12-hour  circular  orbit  with  the  following  instrument  precision: 


Note  that,  even  with  a  four-  to  six-fold  decrease  in  observation 
precision,  the  star-horizon  and  horizon  sensor  navigation  errors  are 
two  to  seven  times  smaller  than  those  predicted  for  the  star-moon 
sensor.  If  the  horizon  sensors  are  an  order  of  magnitude  worse 
(lo  -  30  arcsec),  the  star-horizon  and  horizon  sensor  navigation 
errors  grow  to  350  m  and  230  m,  respectively. 

2.4.2  Star-Horizon  Performance  vs.  Sensor  Orientation 

To  determine  the  optimum  sensor  configuration  for  two  star 
sensors,  the  sensor  bore-sight  direction  was  first  moved  in  azimuth 
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in  the  local  horizon  plane  and  then  in  elevation  for  a  90-minute  and 
a  12-hour  orbit.  Figures  2.10  and  2.11  show  the  effect  of  azimuth 
in  the  90-minute  and  12-hour  orbits,  while  Figures  2.12  and  2.13 
examine  the  effect  of  elevation  in  the  two  cases. 

From  the  figures  it  is  seen  that.,  for  high-altitude 
orbits,  sensors  located  between  45°  and  50°  azimuth  and  as  low  an 
elevation  as  possible  give  the  best  results.  The  lowest  elevation 
possible  coincides  with  the  horizon  sensors  themselves;  thus,  the 
horizon  sensors  alone  produce  the  best  navigation  accuracy  in  the 
absence  of  attitude  errors.  While  the  90-ninute  orbit  was  fairly 
insensitive  to  azimuth  changes,  it  seems  to  be  more  sensitive  to 
elevation  angle  than  does  the  12-hour  orbit.  Interestingly,  a  very 
high  elevation  star  sensor  gives  improved  performance  in  low  orbit, 
and  in  fact,  the  test  program  that  led  to  PADS  had  a  star  sensor 
pair  mounted  60°  from  the  local  horizon,  one  at  0°  azimuth  and  the 
other  at  180°  [Honeywell,  1973]. 

Since  performance  over  the  period  range  90  min  to  5  days 
is  of  interest  to  this  study,  the  star  sensor  configuration  selected 
for  comparison  with  the  other  optical  measurements  in  this  analysis 
was  one  with  star  sensors  located  on  the  spacecraft  in  the  local 
horizon  plane,  ±  1)5°  from  the  T  unit  vector. 

2.5  Optical  Sensor  Performance  vs.  Circular  Orbit  Period 

Fig.  2.1J1  through  2.15  show  the  performance  of  the  optical 
sensors  as  a  functior  of  altitude  or  period.  While  the  star-horizon 
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sensor  shows  an  error  growth  with  altitude  due  to  the  changing 

geometry,  the  space  sextant  exhibits  an  almost  constant  error.  This 
is  due  to  the  distance  of  the  moon  being  so  much  greater  than  the 
satellite  radius.  The  horizon  sensor  shows  the  interesting  feature 
of  a  sharp  drop  in  position  uncertainty  followed  by  nearly  constant 
performance  as  the  altitude  increases.  For  a  constant  sensor  error, 
the  navigation  performance  would  be  expected  to  degrade  with 

increasing  altitude,  as  for  the  star-horizon  sensor.  However,  the 
horizon  sensor’s  apparent  precision  due  to  a  constant  150  meter 
tangent  height  uncertainty  (oh),  falls  rapidly  at  low-to-medium 
altitudes  and  then  drops  more  slowly  as  altitude  increases  (Fig. 
2.16).  Note  that  the  increase  in  apparent  sensor  precision 

parallels  the  navigation  performance  seen  in  Fig.  2. 1 4b ,  i.e.,  as 

altitude  increases  from  near-earth  orbits,  the  rapid  improvement  in 
apparent  sensor  precision  leads  directly  to  corresponding 
improvements  in  navigation  accuracy.  These  results  show  that  if  the 
earth  horizon  can  be  tracked  to  oh  -  150  m,  two  horizon  sensors  can 
provide  much  better  navigation  accuracies  than  the  space  sextant. 
If  the  horizon  detection  is  accurate  to  only  1500  m,  the  resulting 
navigation  error  increases  to  234-350  m  for  a  12-hour  orbit, 
depending  upon  the  mode  used  (horizon  only  or  star-horizon). 

2.6  Sensor  Performance  on  an  Elliptical  Orbit 

Figure  2.17  shows  the  performance  of  the  three  measurement 
types  applied  to  a  12-hour  elliptical  orbit  (e  «  .75)  with  sensors 
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at  ±  ^5°  azimuth  and  0°  elevation.  The  horizon  sensor  must  be 
gimbaled,  however,  and  the  gimbal  mechanism  would  introduce  some 
error  in  the  horizon  definition. 

Comparing  these  results  to  Figure  2.9,  it  is  seen  that  the 
relative  order  of  sensor  performance  stays  the  same  when  applied  to 
an  elliptic  orbit,  but  the  error  magnitude  grows,  especially  for  the 
star-horizon  sensor  pair.  This  error  growth  could  be  reduced  by 
selecting  a  lower  elevation  sensor  direction  as  altitude  is 
increased,  but  that  would  lead  to  a  much  more  complicated  system 
unless  the  star  and  horizon  sensor  pair  were  mounted  on  the  same 
gimbaled  platform.  If  that  were  the  case,  the  relative  orientation 
of  the  two  sensors  would  remain  the  same,  and  gimbal  inaccuracies 
could  probably  be  reduced  considerably. 

2.7  Comparative  Sensor  Performance  for  Circular  Orbits 

The  results  of  this  analysis  are  summarized  in  Table  2.1, 
where  the  matrix  sensors  are  operated  in  either  the  star-horizon 
angle  mode  or  as  horizon  sensors  alone.  The  navigation  accuracies 
are  representative  of  orbits  above  300-minute  period. 
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TABLE  2.1 

PREDIC  iD  NAVIGATION  ERROR  DUE  TO  SENSOR  ERROR 


Sensor 

Sensor  Precision 

Navigation  Error 

[arcsec  J 

[meters  j 

Space  Sextant 

.5 

150 

Matrix  Sensors 

Star-Horizon 

2 

40-180 

1 

19-90 

.5 

9-45 

Horizon 

2 

20 

1 

11 

.5 

6 

It  should  be  emphasized  that  only  sensor  error  was 
considered  and  that  the  fixed  sensor's  results  would  be  degraded  by 
attitude  errors.  This  is  especially  true  for  the  horizon  sensor, 
since  the  vehicle-determined  local  vertical  forms  the  measurement 
reference.  Note  that  the  optical  sensors  may  also  form  the  attitude 
determination  sensor  system,  as  indicated  in  the  proposed 
utilization  of  the  space  sextant.  Correlation  between  attitude  and 
navigation  solutions  may  be  a  problem,  but  since  the  attitude  filter 
requires  high  frequency  updates  while  the  navigation  filter  may 
require  less  frequent  update,  the  two  requirements  may  be  handled  as 
two  separate  problems  as  long  as  the  navigation  system  is  aware  of 
the  time  history  of  the  attitude  solution. 

2.8  Refraction  Errors 

The  largest  error  sources  inherent  in  the  star-earth 
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horizon  measurement  of  star  refraction  are  the  sensor  itself  and  the 
atmospheric  modeling  required  to  predict  refraction  angles.  In 
studies  recently  completed  for  the  Naval  Surface  Weapons  Center 
(NSWC)  by  Clynch  [l979  and  1981b],  a  ray  tracing  program  was  used  to 
compute  refraction  of  a  star  ray  passing  through  earth's  atmosphere. 
The  model  for  refract! vity  was  from  Owens  [ 1 967  ]  where  the 
refract! vity ,  N,  is  a  function  of  the  refractive  index, 

N  5  I06(n-l)  (2.25) 


but  the  refractive  index,  n,  of  a  gas  must  be  determined  by  the 
molecular  density  and  relative  polarizability.  This  refractive 
index  is  obtained  from  the  Lorentz-Lorenz  equation. 


where  Ri  and  are  the  specific  refractivity  and  density  of  the 
ith  component.  Assuming  that  dry,  C02~free  air  can  be  treated  as 
single  component,  equation  (2.26)  can  be  approximated  by 

-  R  P  (2.27) 

xf+2 

From  (2.27),  it  is  seen  that 


n 


2 


1  +  2  R  p 

4  ^  A 

1  "  R  p 


(2.28) 
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Thus, 


(2.29) 


Note  that  R  in  Eq.  (2.29)  is  a  constant  for  a  homogenous 
dry  air  mixture,  but  that  density,  p,  is  not  a  constant  nor  is  it  a 
well-behaved  function  of  altitude.  Density  is  usually  modeled  as  an 
exponential  function  of  altitude, 


(2.30) 


where 


P0  -  bottom  altitude  density 
a  -  altitude 
Hg  *  scale  height 

but  this  exponential  model  breaks  down  when  attempting  to  model 
large  altitude  bands.  To  extend  the  range  of  the  equation,  the 

altitude  is  broken  into  several  altitude  bands  of  constant 
thickness.  The  constants  in  Eqs,  (2.29)  and  (2.30),  R  ,  pQ  and  H, 
also  vary  as  a  function  of  location  and  of  local  weather  conditions; 
thus,  location  errors  and  delays  result  in  errors  in  computed 

refraction  index  and  thus  in  the  computed  refraction  angle  for  a 
given  ray  altitude.  The  computed  refraction  angle,  6,  in  a  slab  is 
given  by  a  continuous  form  of  Snell's  law  in  the  case  where  a  slab 
is  defined  as  a  region  in  which  n  has  a  small,  constant  gradient 
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[Kelso,  1964] 

n  sin  8  -  constant  -  nQ  sin  eQ  (2.31 ) 

The  computer  program  used  by  Clynch  propagates  the 

starlight  through  the  atmosphere  and,  by  constructing  a  series  of 

*. 

slabs  of  constant  refractive  index  gradient,  computes  the  bending  of 
the  ray  as  it  moves  from  one  slab  through  another.  In  this 
analysis,  the  log  of  density  vs.  altitude  data  obtained  from  the 
National  Oceanic  and  Atmospheric  Administration  were  fit  with  cubic 
polynomials  with  an  average  relative  error  of  0.5  percent.  When 
these  models  were  then  used  to  predict  atmospheric  density  and  then 
RTA  (Fig.  2.7),  it  was  found  that  the  error  in  RTA  was  approximately 
proportional  to  the  relative  error  in  density  at  the  ray  tangent 
point  with  a  1  percent  density  error  producing  a  100  m  RTA  error. 
When  data  one  to  four  days  apart  were  processed,  the  results  showed 
the  effect  of  data  aging.  Table  2.2  li3ts  typical  errors  of  Ap/p 
for  several  one-  and  four-day  spans  in  1979.  These  values  are 
computed  near  the  north  pole,  and  latitudes  below  which  density 
changes  of  si  5  percent  were  observed  are  indicated.  Winter  pole 
density  errors  are  very  large  compared  to  summer,  but  the  error 
quickly  drops  as  warmer  weather  approaches.  Even  in  the  worst 
month,  however,  5  percent  error  or  less  is  observed  over 
approximately  60  percent  of  the  earth's  surface. 
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TABLE  2.2 


EFFECT  OF 

ATMOSPHERIC 

DATA  AGEING  AT 

20  KM  ALTITUDE 

Date 

Age  (Days) 

Ap/p(*) 

Latitude 

9  -Ian 

1 

7.5 

40° 

12  Jan 

4 

27.5 

40° 

21  Jun 

1* 

4.4 

90  0 

Clynch  states  that  when  data  are  aged  for  five  days  or  more,  they 
are  comparable  to  climatology  predictions;  thus,  density  predictions 
covering  the  quiet  60  percent  of  the  earth  should  be  accurate  to 
approximately  5  percent. 

These  errors  are  those  expected  from  a  complex  ray  tracing 
computer  program  operating  with  aged  atmospheric  data.  In  the 
onboard  navigation  scenario,  it  is  not  possible  to  implement  this 
ray  tracing  algorithm,  and  a  reduced  order  model  mu3t  be  used  that 
includes  an  atmosphere  generated  onboard  with  a  minimum  of  updates. 
Clynch  [ 1981  a]  proposes  an  approximation  to  the  ray  path  that 
matches  the  actual  density  and  density  gradient  only  at  the  ray 
tangent  point  (rtp).  Altitude  (a)  near  the  RTP  is  approximated  as  a 
function  of  distance  from  the  RTP,  (s),  by 


(2.32) 


where 
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rr  -  Re  +  RTA 


The  refraction  angle  0  is  then  expressed  as  the  integral  over  the 
path  length 


e(RTA)  -  j 


(2.33) 


The  actual  path  length  derivative  de/ds  is  a  function  of  conditions 
in  each  slab,  but  when  the  actual  optical  path  is  replaced  by  a 

straight  line  through  the  RTP  and  tangent  to  the  ray  there,  the 

2 

integral  becomes  a  Gaussian  integral  (containing  e-x  dx  terms). 


e(RTA)  -  3fe(l).P. 
2H 


r  K-) 

L  \2rrHs/ 


3R(l )rh 


(2-3«) 


where  R ( A )  is  the  specific  refractivity  of  dry  air  at  the 
wavelength,  A  ,  of  interest,  and  p  is  evaluated  at  the  RTP. 

If,  in  addition  to  the  assumption  that  curvature  of  the 
path  takes  place  very  near  the  RTP,  it  can  also  be  assumed  that 
conditions  about  the  RTP  are  stable,  then  data  on  the  atmosphere  at 
the  geometric  tangent  point  (the  geometric  intersection  of  tangent 
lines  from  the  star  and  spacecraft  •  GTP)  can  be  used  to  determine 
H  and  0(rta).  RTA  is  then  approximated  by  geometric  tangent 
altitude  (OTA ) .  This  greatly  simplifies  the  problem  since  GTA  is 
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determined  by  star ,  earth  and  spacecraft  positions  alone.  The 
equation  relating  e(GTA)  and  9(rta)  is 

9(RTAj  -  0(GTA)  exp  [-B0[RTA)]  (2.35) 

where 


Sp  -  distance  from  observer  to  point  of  ray  tangency 

Since  Eq.  (2.35)  is  transcendental  in  0(rta),  large  differences  in 
atmospheric  conditions  (hs  and  p)  between  GTA  and  RTA  would 
require  an  iterative  solution. 

Testing  this  simplified  model  for  o(rta)  vs.  the  complex 
atmosphere  in  the  ray  tracing  program  showed  a  relative  error  of 
approximately  2  percent  in  the  GTA  range  25  to  iJ5  km.  This 
corresponds  to  an  error  of  2.JJ  to  .16  arcsec  in  o(rta)  .  These 
results  indicate  that  the  contribution  of  atmospheric  modeling  to 
the  CCD  sensor  error  budget  is  on  the  order  of  2  arcsec  plus  the 
error  due  to  not  having  current  atmospheric  data.  If  monthly 
atmospheric  data  are  transmitted  to  the  vehicle  and  if  that  data  can 
then  be  extrapolated  to  produce  less  than  an  8  percent  density 
error,  then  the  refraction  model  would  contribute  approximately 
1000  m  to  the  error  in  GTA. 

Even  a  total  error  in  density  of  15  percent,  giving  a  GTA 
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error  of  1500  m  would  result  in  a  12-hour  orbit  navigation  error  of 
only  350  m  with  two  star-horizon  CCD  sensor  pairs.  Lower  orbits 
would  feel  that  GT£  er^or  more  acutely,  and  higher  orbit  accuracy 
would  be  better.  Before  realistic  error  budgets  can  be  produced, 
however,  more  work  must  be  done  to  quantify  errors  in  long-term 
density  prediction. 
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3. 1  Introduction 

The  results  of  Sections  2.7  and  2.8  show  that  the  optical 
sensors  being  developed  will  produce  navigation  accuracies  on  the 
order  of  300  meters.  Several  satellite  programs,  however,  require 
errors'  of  less  than  50  m  and  are  therefore  potential  GPS  users  if 
GPS  can  be  made  less  dependent  on  ground  tracking  and  control.  This 
is  the  purpose  of  the  remainder  of  this  report  —  to  investigate  the 
feasibility  of  autonomous  navigation  of  GPS  vehicles  using 
satellite-to-satellite  range  and  integrated  doppler  Information  in  a 
reasonably  sized  onboard  navigation  processor.  The  approach  taken 
is  to  first  determine  constellation  selection  effects,  then  to 
analyze  filter  model  requirements  with  consider  covariance 
techniques  and  to  propose  an  onboard  filter  design  and  evaluate  its 
error  propagation  characteristics. 

3*2  GPS  Constellations 

The  Global  Positioning  System  is  composed  of  a  spaceborne 
segment  consisting  of  a  variable  number  of  satellites  in  12-sidereai 
hour  orbits,  a  control  segment  consisting  of  a  Master  Control 
Station  (MCS } ,  and  four  remote  Monitor  Stations  (ms)  and  a  user 
segment  consisting  of  many  different  user  receiver  sets  built  by 
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several  different  manufacturers.  The  number  of  satellites  to  be 
placed  in  the  constellation  varies  with  congressional  budget 
decisions.  Descriptions  of  the  system  and  its  several  facets  abound 
in  the  literature,  and  an  excellent  collection  of  papers  can  be 
found  published  together  by  the  Institute  of  Navigation. 

The  number  of  satellites  currently  planned  is  18,  and  two 
different  constellations  have  been  analyzed  by  the  Joint  Program 
Office  (JPO)  [Book,  et  al . ,  1 980 ] .  The  first  of  these  is  a  modified 
three-plane  2iJ-satellite  configuration  in  which  two  of  the  eight 
satellites  in  each  plane  are  eliminated  leaving  an  unsymmetric  18- 
satellite  constellation.  This  is  termed  the  3x6  configuration 
since  three  planes  of  six  satellites  each  are  employed.  The  second 
candidate  constellation  consists  of  six  planes  with  three  satellites 
each  in  a  symmetric  pattern,  i.e.,  a  6  x  3  constellation.  The  orbit 
elements  for  these  constellations  are  given  in  Table  3.1 • 

According  to  Book,  et  al.,  the  unsymmetric  3x6 
constellation  has  a  geometric  performance  with  respect  to  the  ground 
almost  as  good  as  the  symmetric  6x3  configuration,  where  GPS 
geometric  performance  is  determined  by  examining  the  trace  of  a  unit 
covariance  matrix  of  the  user  position  error  as  follows  [Milliken 
and  Zoller,  1 980 ] . 

If  the  user  receives  information  from  at  least  four  GPS 
satellites,  he  can  estimate  his  position  and  his  clock  error  and  the 
error  in  these  estimates.  He  measures  a  pseudo-range  to  each 
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satellite: 


-  Pi  +  cAtAi  +  c(Atu  -  At  )  +  ci 


(3.1) 


where 

*•  Pi  »  pseudo-range  to  satellite  i 
Pi  -  geometric  range 

At  *  satellite  i  clock  offset  from  a  common  time  scale 

Atu  -  user  clock  offset 

At.  .  propagation  delays,  etc. 

Ai 

■  random  measurement  noise 


If  AtAi  is  determined  or  adequately  modeled  and  Atg  known  for 

each  satellite,  then  four  Pj  measurements  to  four  GPS  vehicles 
provide  information  to  solve  for  user  position  (x,y,z)  and  clock 


bias  (cAtu).  Since  pi 


1/2 


(*„  -x)2  +  (y„  -y)2  +  (z_  -z )‘ 


bi 


,  the 


measurement  partial  derivatives  are  given  by 
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(3.2) 


Now,  assuming  a  unit  error  variance  in  ,  a  user  information  matrix 


can  be  computed  from 


"2  •••  »I]  [I] 

and  the  corresponding  unit  covariance  matrix  is 

C0Vu  -  inf'1.  (3.^) 

This  is  termed  a  unit  covariance  here  because  the  unit  range 
variances  do  not  reflect  reality  but  produce  user  position  sigmas 
that  can  then  be  multiplied  by  actual  system  errors  to  give 
realistic  user  errors.  The  trace  of  the  COVu  matrix  describes  the 
variance  of  the  user's  position  error  and  clock  bias  given  a  unit 
variance  in  each  range  measurement.  Tr[c0Vu]  iS  thus  called  the 
Geometric  Dilution  of  Precision  (GDOP)  since  it  describes  a  four¬ 
dimensional  geometri.  error.  Other  similar  definitions  follow: 

PDOP  -  [C0Vu(l , 1 )  +  C0VU(2,2)  +  C0VU(3, 3) ]  (3-5) 

=  Position  Dilution  of  Precision 

HDOP  -  [COV  (2,2)  +  COV..  (3,3)]  (3.6) 

RTN  URTN 

*  Horizontal  Dilution  of  Precision 


(3.3) 


k  i  4 
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VDOP  -  COV  (l,l)  (3.7) 

RTN 

B  Vertical  Dilution  of  Precision 
TDOP  -  COVu(l|,U)  (3.8) 

*  Time  Dilution  of  Precision 

Note  that  GDOP  requires  the  inversion  of  a  4  x  4  symmetric 

matrix,  and  PDOP  requires  the  inversion  of  a  3  x  3  matrix.  It  will 

be  seen  in  GDOP  and  PDOP  plots  that  they  provide  essentially  the 

same  information;  thus,  a  user  can  save  processor  time  by  using  PDOP 

in  selecting  the  optimum  set  of  four  GPS  satellites.  It  would  be 

even  more  beneficial  if  the  trace  of  the  user  information  matrix 

instead  of  COVu  could  he  used  to  determine  optimum  GPS  satellite 

selection.  But,  as  noticed  by  Fang  [l980],  each  row^^  partial  in  the 

H  matrix  is  a  unit  vector  pointing  to  GPS^  while  the  fourth 

T 

component  is  1.  When  the  product  H  H  is  formed,  the  resulting 
trace  for  n  satellites  can  be  obtained  through  the  commutative 
property  of  the  trace, 

Tr(HTH)  -  7r(HHT)  (3-9) 

but, 

TrtHHT)  •  j,  {S  *  (-it)  *  iS  * '  • 2n  (3-’o) 
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Therefore,  the  trace  of  the  user  information  matrix  is  a  constant 
equal  to  2n  for  n  satellites  and  gives  no  information  about  the 
geometric  attributes  of  the  configuration. 

3 • 3  Constellation  Comparison 

To  determine  the  effect  of  constellation  selection  (3x6 
or  6x3)  on  GPS  vehicle  navigation,  a  relative  geometry  computer 
program  was  written  to  step  through  one  12-hour  GPS  cycle.  Each  10 
minutes,  satellite  positions  were  computed  and  used  to  determine 
GDOP  and  PDOP  for  any  desired  GPS  user  spacecraft.  In  addition, 
satellite  visibility  times  were  accumulated  so  a  total  percentage  of 
the  orbit  visible  to  each  satellite  was  available.  Figure  3*1  shows 
GDOP  and  PDOP,  along  with  number  of  satellites  visible  to  the  GPS  #1 
spacecraft  for  the  symmetric  6x3  constellation.  Elevation  limits 
of  -28°  to  -76°,  as  recommended  by  Chao  [ 1 981 ] ,  were  used  in  the 
visibility  check,  and  all  satellites  visible  were  included  in  the 
calculations.  It  was  decided  to  use  all  visible  satellites  in  the 
geometry  calculate  3  because  the  GPS  crosslink  (l^  at  1381  MHz) 
will  allow  each  satellite  to  transmit  for  approximately  1.5  seconds 
every  36  seconds  [Barr,  1981 ]  and  the  navigation  algorithm  would 
benefit  from  clock  and  position  information  from  all  available 
vehicles  rather  than  selecting  an  optimum  set  of  H.  GDOP  and  PDOP 
follow  similar  but  offset  paths  for  other  vehicles,  but  the  GDOP 
patterns  are  identical  for  any  two  vehicles  when  corrected  for 


DGP  &  «SflTS  VISIBLE/2 


Figure  3.1.  GPS  Geometric  Performance,  6x3  Symmetric 
Constellation,  Satellite  No.  1 
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satellite  phase  difference.  Note  that  there  is  an  inverse 
correlation  between  GDOP  and  number  of  satellites  visible,  with  the 
minimum  GDOP  (2.5}  corresponding  to  14  vehicles  visible  and  the 
maximum  ("4.2)  corresponding  to  11  visible. 

Figure  3.2  shows  the  same  information  for  the  unsymmetric 
3x6  constellation,  satellite  #1 .  This  constellation  does  not 
produce  similar  GDOP  patterns  for  different  vehicles,  and  the 
minimum  values  (“2.5  and  “2.8)  and  maximum  values  (“4.7  and  4.o)  are 
also  inconsistent  among  users.  Table  3.2  shows  the  minimum  and 
maximum  values  over  one  orbit  for  the  18  satellites.  Note  that  most 
of  the  values  are  higher  than  those  encountered  in  the  6x3 

constellation  and  that  all  of  the  average  values  are  higher  in  the 
3x6  case. 

Figure  3.3  shows  the  percent  of  the  orbit  during  which 
each  GPS  se'ellite  is  visible  to  vehicle  # 1  in  the  6x3 

constellation.  The  values  range  from  a  low  of  about  38?  to  several 
100$  cases.  Figure  3.4  depicts  the  same  information  for  vehicle  //I 

in  the  3x6  constellation.  The  main  difference  is  that  some 

satellites  in  the  3x6  constellation  are  invisible  to  each  other 
(satellite  # 1  cannot  see  #2  and  #5),  potentially  weakening  the 
position  and  clock  estimatloon  results  of  the  whole  system.  For  this 
reason  and  also  because  of  better  GDOP  performance,  the  symmetric 
6x3  constellation  is  desirable  from  a  satellite-to-satellite 


tracking  viewpoint. 
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Figure  3.2.  GPS  Geometric  Performance,  3x6  Unsymroetric 
Constellation,  Satellite  No.  1 


TABLE  3.1.  GPS  ORBIT  ELEMENTS 


A.  Unsymmetric  3x6  Constellation 
Semi-major  axis:  26575.6  km 
Eccentricity:  0.0 

Inclination:  55  deg 

Relative  Ascending  Node:  0,  120,  240  deg 

Argument  of  Perigee:  0  deg 

Mean  Anomaly  at  Epoch: 

Plane  1:  0,  45,  90,  135,  270,  315  deg 

Plane  2:  20,  65,  110,  155,  200,  245  deg 

Plane  3:  40,  175,  220,  265,  310,  355  deg 


B.  Symmetric  6x3  Constellation 
Semi-major  Axis:  26575.6  km 

Eccentricity:  0.0 

Inclination:  55  deg 

Relative  Ascending  Node:  0,  60,  120,  180,  2 40,  300  deg 

Argument  of  Perigee:  0  deg 

Mean  anomaly  jf  Epoch: 

Plane  1 :  0,  120,  240  deg 

Plane  2:  ;,0,  i60,  280  deg 

Plane  3:  80,  200,  320  deg 

Plane  4:  120,  240,  0  deg 

Plane  5:  160,  280,  40  deg 

Plane  6:  200,  320,  80  deg 


3.4  GPS  Error  Sources 

The  consider  covariance  analysis  that  follows  in  Chapter  4 
provides  realistic  satellite  position  errors  when  given  realistic 


TABLE  3.2 

MINIMUM,  MAXIMUM  AND  AVERAGE  GDOP  VALUES 

A.  Unsymmetric  3x6  Constellation 

Satellite:  1234567 

8  9 

Min  GDOP: 

2.37 

2.47 

2.45 

2.36 

2.46 

2.58 

2.38 

2.46 

2.60 

Max  GDOP: 

4.40 

4.04 

4.72 

4.29 

4.21 

4.70 

4.29 

4.83 

4.51 

Avg  GDOP: 

3-30 

3.25 

3.48 

3-19 

3.24 

3.46 

3.30 

3.62 

3.51 

Sa.  llite: 

10 

11 

12 

13 

14 

15 

16 

17 

18 

Min  GDOP:  2.65  2.66  2.29  2.58  2.29  2.59  2.38  2.38  2.64 

Max  GDOP:  4.54  4.62  4.20  4.23  10  4.14  4.48  4.58  4.83 

Avg  GDOP:  3.44  3-60  3.28  3.23  3.19  3-37  3-25  3-35  3.49 


B.  Symmetric  6x3  Constellation 


Satellite: 

All 

Min  GDOP: 

2.37 

Max  GDOP: 

4,25 

Avg  GDOP: 


2.88 
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values  of  the  1o  errors  expected  to  be  seen  by  the  satellite.  The 
purpose  of  this  section  is  to  describe  the  error  sources  in  the  GPS 
autonomous  navigation  scenario  and,  where  necessary,  to  determine 
values  or  bounds  on  the  values  expected  to  be  experienced  by  the 
satellite.  The  results  of  the  consider  analysis  will  then  indicate 
which  error  source  effects  must  be  modeled  or  estimated  in  the 
onboard  filter  algorithm.  Some  of  these  effects  are  unobservable  by 
satellite-to-satellite  tracking,  so  these  error  sources,  if  the 
effects  are  large  enough,  will  be  added  to  the  onboard  algorithm 
without  considering  them  in  the  covariance  analysis. 

The  GPS  autonomous  navigation  problem  contains  several 
error  sources  in  common  with  a  user  on  the  ground  and  some  unique 
errors  of  its  own.  A  ground  user  experiences  the  following  error 
sources  [General  Dynamics,  1978]: 

a.  user  clock  bias  and  drift, 

b.  satellite  clock  bias  and  drift, 

c.  receiver  movement  during  signal  transit, 

d.  satellite  ephemeris, 

e.  relativistic  effects, 

f.  antennae  offsets, 

g.  receiver  signal  delay, 

h.  time  tagging, 

i.  ionospheric  delay, 

j.  tropospheric  delay. 
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A  GPS  satellite  receiver  located  above  the  ionosphere  sees  errors 
due  to 

a.  user  clock  bias  and  drift, 

b.  satellite  clock  bias  and  drift, 

c. -  receiver  movement  during  signal  transit, 

d.  satellite  ephemeris, 

e.  • relativistic  effects , 

f.  antennae  offsets, 

g.  receiver  signal  delay, 

h.  time  tagging, 

i.  earth  geopotential , 

j.  n-body  gravity, 

k.  solar  radiation  pressure, 

l.  vehicle  thrusting  or  outgassing, 

m.  earth  polar  motion  and  angular  velocity. 

In  this  study,  it  is  assumed  that  pre-processing  of 
pseudo-range  and  integrated  doppler  measurements  is  performed  and 
that  residual  errors  due  to  receiver  movement,  relativity,  antennae 
offsets,  receiver  signal  delay  and  time  tagging  can  be  modeled  or 
corrected  for  with  the  remaining  measurement  errors  expressed  as 
uncorrelated,  zero-mean,  Gaussian  errors.  Simulation  of  the 
remaining  errors  is  described  in  the  following  sections. 


73 


3.^.1  GPS  Clock  Errors 

Operational  GPS  satellites  are  planned  to  contain  two 

rubidium  and  two  cesium  beam  frequency  standards  with  the  possible 

addition  of  a  hydrogen  maser  [Payne,  1982].  To  characterize  the- 

errors  inherent  in  these  precise  time  standards,  it  is  convenient  to 

measure  clock  stability  in  terms  of  the  Allan  variance  of  its 

fractional  frequency  error,  ■*—  .  The  error  in  the  GPS  clock  after 

o 

some  elapsed  time  from  update  is 


T(f)  -  T(fo) 


where 


fQ  -  2mv0  =  nominal  frequency  (rad/sec] 
f  »  time-varying  true  frequency 
Af  »  frequency  offset 
f  -  frequency  drift 

6f(t)  »  time-varying  random  frequency  error 


(3.11) 


Differentiating  yields 


o 


o 


o 


o 


(3.12) 
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-20  db/decade  slope  versus  frequency.  In  addition,  at  very  low 
frequencies  (long  sample  times,  t ) ,  the  integral  of  white  noise 
causes  a  -*40  db/decade  slope.  A  typical  Allan  variance  curve 
showing  GPS  specifications  is  given  in  Figure  3.5,  along  with  actual 
on-orbit  clock  values  [Payne,  1982]. 

3. 4.1.1  GPS  Clock  Error  Simulation 

Since  Allan  variance  is  widely  used  as  a  measure  of  atomic 
frequency  standard  stability,  time  error  simulation  using  an  input 
value  for  Allan  variance  is  a  necessity  for  GPS  error  analyses. 

To  generate  a  frequency  error  signal,  y(t)  ,  and  its 
integral,  the  phase  error,  the  Allan  variance  curve  is  described  in 
the  time  domain  by  the  power  series 

°y(T)  -  l  K  TS  (3.17) 

8  e 

where  8  *  -2,  -1,  0,  +1,  +2  and  each  value  of  8  dominates  in  a 
region  to  Tj  .  In  the  frequency  domain,  a  similar  series  can 

be  used  to  compute  the  one-sided  spectral  density, 

Sy(f)  ‘  l  hQ  fa  (3.18) 

a 

with  a  having  the  same  values  as  8  .  Similar  to  6  ,  each  value 
of  a  dominates  over  a  certain  frequency  range.  The  constants, 
and  h^  define  the  level  of  the  time  or  frequency  clock  error 
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segment. 

To  simulate  y(-r)  ,  the  Allan  variance  levels  and  corner 
frequencies  are  specified.  These  values  determine  the  Bode  plot  of 
the  one-sided  spectral  density.  From  the  Bode  plot,  spectral 
shaping  produces  a  dynamic  system  function  [Laning  and  Battin, 
1977].  An  inverse  Laplace  transform  then  generates  time  domain 
linear  differential  equations.  The  solution  of  these  equations 
produces  y(-r)  and  its  integral. 

For  white  frequency  noise  and  integrated  white  frequency 
noise,  the  inverse  Laplace  transform  is  straightforward.  The  power 
spectral  density  function  corresponding  to  integrated  white  noise 
(frequencies  below  f^  or  averaging  times  above  T2)  is 


sy(s) 


b2(s2  +  a2o>2) 
a2(s2  +  ui2) 


(3.19) 


where  3  is  the  complex  frequency  0  +  jio  and  m  »  f. 


This  is  the  spectral  density  of  a  Gaussian  white  noise  process 
driving  a  linear  system.  It  can  then  be  expressed  as 

Sy(s)  -  h(s)  H* (s )  Su(s )  (3.20) 

where  h(s)  is  the  linear  system  transfer  function,  H*(s)  is  its 
complex  conjugate  and  Su(s)  is  the  spectral  density  of  the  white 
noise  process.  Since  Su(s)  is  constant,  assume  it  to  be  unity, 
implying  unit  variance,  then 
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h(s) 


U  (s  +  ao>1 ) 

1  (a  ♦  «,') 


and  the  transfer  function  for  the  output  signal  is 


Y(s)  -  H(s)  u(s)  -  H(s ) . 


(3.21) 


(3.22) 


Taking  an  inverse  Laplace  transform  produces  the  linear  differential 
equation 

x(t)  «  -0)^(7)  +  (a-l)oj1u(x)  (3.23) 

then 

y(t)  -  j  [x(t)  +  u(x)J  (3.2H) 

where  u(t)  is  the  output  of  a  unit  variance  Gaussian  random  number 
generator . 

Flicker  noise,  however,  cannot  be  generated  by  one  inverse 
Laplace  transform  since  the  spectral  density  curve  corresponds  to  a 
transfer  function  of  the  form 


H(s ) 


(s)1/2 


1 1I/2 


(3.25) 


for  which  a  finite-order  state  representation  cannot  be  constructed 
which  will  generate  the  system  output  [Meditch,  1975).  Several 
approaches  to  approximate  y(t)  for  flicker  noise  have  been 
proposed,  and  Meditch  describes  a  method  by  Barnes  and  Jarvis  [ 1 97 1  ] 
that  efficiently  models  1/(s)1/2  by  a  cascade  of  lag  networks. 
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The  technique  consists  of  approximating  the  desired  transfer 
function  with  slope  -10  db/decade  with  n  stages,  each  consisting 
of  a  -20  db/decade  section  and  a  white  noise  (constant  spectral 
density)  section.  For  simulation  purposes,  Meditch  states  that  a 
choice  of  n  -  4  gives  a  reasonable  and  efficient  approximation  to 
the 'desired  shape. 

Once  the  series  approximation  to  1/(s)^2  is  complete, 
the  inverse  Laplace  transform  gives  a  linear  differential  equation 
for  each  stage.  These  equations  then  form  an  n-vector  linear 
differential  equation  whose  solution  provides  x(t)  such  that  y(x) 
for  averaging  times  between  and  can  be  simulated. 

As  seen  in  Figure  3-5,  GPS  clocks  are  exhibiting  Allan 
variances  of  lO"1^  for  t  -  1  hour  to  10-1^  for  r  >  2  days.  These 
values  give  stabilities  of  10  J  to  10  ns/s  and  are  used  in  the  GPS 
navigation  error  model  as  expected  accuracies  of  operational  cesium 
clocks.  It  is  assumed  that  hydrogen  masers  are  one  order  of 
magnitude  better  [Kartaschoff ,  1978,  p.  62].  Note  that,  from  Figure 
3.5,  another  order  of  magnitude  improvement  in  6f/f  would  allow 
the  clocks  to  run  independently  for  140  days  with  a  user  equivalent 
range  error  (uERE)  of  approximately  90  meters. 

3.4.2  Satellite  Ephemerls  Errors 

GPS  ephemeris  errors  in  the  ground  tracking  mode  are 
required  to  be  on  the  order  of  1-10  m.  Current  estimates  of 
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position  errors  indicate  that  12-hour  periodic  errors  of  2,  10  and 
6m  in  the  RTN  directions  are  combined  with  a  secular  tangential 
error  growth  of  3  m/day  [Anderle,  1 980 ] .  This  short-term  linear 
growth  becomes  quadratic  in  the  long  term.  GPS  satellite-to- 
satellite  tracking  errors,  however,  have  not  been  quantified  for 
long-term  operation.  This  is  one  goal  of  this  study.  Liu's  [ 1 981 ] 
analysis  of  satellite  ephemeris  errors  over  two  weeks  of  GPS-GPS 
tracking  shows  an  80  m  secular  growth  and  a  six-hour  20  m  periodic 
term  due  to  solar  radiation  pressure  errors. 

In  this  study,  one  satellite  at  a  time  is  analyzed  and  it 
is  assumed  that  initial  RTN  error  sigmas  are  2,  10  and  6  m, 
respectively,  with  a  secular  tangential  growth  of  3  m/day  for  the 
other  vehicles.  These  values  are  used  in  the  consider  analysis  to 
determine  navigation  errors  of  satellite  //I. 

3. ^-3  Earth  Geopotential  Errors 

One  of  the  goals  of  this  study  is  to  determine  the 
accuracy  and  size  (order  and  degree]  of  the  onboard  geopotential 
required  for  accurate  GPS  navigation.  In  addition  to  being  subject 
to  secular  perturbations  due  to  J2,  the  GPS  12-hour  orbit  is 
resonant  with  the  harmonic  coefficients  of  degree  2  and  *1,  so  errors 
in  their  values  and  errors  due  to  their  absence  in  the  onboard 
filter  will  be  magnified  over  any  long  prediction  interval. 

Table  3.3  from  Wagner  and  Lerch  [ 1 978]  describes  the 
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estimated  error  in  the  GEM  8  earth  model  obtained  by  comparing  the 
earth  model  predictions  with  new  observations  not  included  in  the 
model  formulation.  These  errors  are  then  used  in  the  analysis  of 
geopotential  error  effects  on  GPS  orbits  in  Chapter 


TABLE  3.3 

ESTIMATED  GEM  8  GEOPOTENTIAL  ERRORS  (  x  10"9} 
Degree  ( £) 
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3«  ^ ^  Solar  Radiation  Pressure 


As  stated  in  Section  3-4.2,  solar  radiation  errors 

considered  by  Liu  caused  a  20  m  periodic  satellite  position  error. 

Liu  assumed  a  1o  radiation  pressure  error  of  10J,  for  an 

-10  2 

acceleration  uncertainty  of  10  m/s  .  Recent  discussions  with  the 
NSWC  personnel  who  determine  GPS  reference  orbits  indicate  that 
radiation  pressure  coefficient  errors  (lo)  of  approximately  1>  have 


been  observed. 


These  figures  provide  a  range  for  the 


errors 
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considered  in  this  analysis. 

3. ^.5  Vehicle  Thrusting  or  Outgassing 

The  GPS  vehicle  is  subject  to  periodic  gas  jet  thrusting 
to  maintain  or  change  orbital  elements  and  to  dump  excess  momentum 
from  the  attitude  control  system  (ACS)  momentum  wheels.  The  orbit 
adjusts  are  planned  and  occur  infrequently,  but  when  gas  jet 
momentum  dumping  occurs,  it  is  performed  automatically  by  the  ACS. 
If  operational  satellites  employ  gas  jet  dumping,  satellite 
navigation  performance  would  suffer  dramatically  each  time  gas  jet 
firings  occur  because  of  thrust  imbalance  and  possible  plume 
impingement.  It  is  expected,  however,  that  magnetic  dumping  of 
momentum  will  be  accomplished  by  the  onboard  ACS  processor.  Tests 
on  current  GPS  vehicles  indicate  that  magnetic  momentum  dumping  is 
successful  [Ferguson  and  Kronke,  1 980 ] ,  so  gas  jet  thrusting  is 
necessary  only  for  orbit  maneuvering. 

Outgassing  is  a  phenomenon  that  does  not  lend  itself  to 
easy  prediction.  Phase  I  GPS  vehicles  on  orbit  appear  to  be 
experiencing  an  acceleration  along  the  spacecraft  solar  panel  axis 
on  the  order  of  +  10-12  to  ±  10-1^  m/s2.  The  cause  of  this 
acceleration  is  possibly  due  to  unmatched  radiators  on  opposite 
sides  of  the  vehicle,  and  this  thrust,  while  not  strictly 
outgassing,  has  a  form  similar  to  that  caused  by  the  boiling  off  of 
volatile  gases.  If  analysis  confirms  that  the  operational  vehicles 
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may  experience  unmodeled  thrusting,  then  the  onboard  filter  must 
include  these  unmodeled  accelerations  in  the  state  vector.  Since 
the  existence  and  form  of  these  accelerations  in  the  operational 
system  are  unknown,  they  are  neglected  in  this  analysis. 

3.4.6  Earth  Polar  Motion  and  Angular  Velocity 

Even  if  the  GPS  satellites  could  navigate  accurately  with 

*  ( 

respect  to  each  other  and  contained  precise  models  of  the  other 
perturbing  forces,  the  constellation  would  still  drift  from  the 
earth-centered,  rotating  frame  (ECR)  in  which  user  positions  are 
defind.  The  largest  error  source  between  an  inertial  earth-centered 
frame  (ECl)  and  the  ECR  frame  is  the  angular  velocity  of  the  earth, 
with  smaller  errors  caused  by  polar  motion.  The  integral  of  angular 
velocity  errors  is  the  time  difference  UT1-UTC.  Currently,  it  has  a 
yearly  drift  of  approximately  1  second,  the  well-known  leap  second 
correction. 

UTC  (Coordinated  Universal  Time)  represents  a  "uniform 
time  scale"  and  is  obtained  by  applying  a  fixed  offset  of 
32.184  seconds  to  an  international  atomic  time  scale  (TAl) 
maintained  by  the  Bureau  International  de  l’Heure  (BIH)  in  Paris 
[The  Astronomical  Almanac,  1983].  The  difference,  UT1-UTC,  plus  the 
record  of  leap  seconds  then  is  a  measure  of  the  changes  in  earth's 
rotation  rate.  UT1-UTC  is  published  for  5-day  intervals  by  the  BIH, 
and  these  tables  provide  the  raw  data  by  wnich  UT1-UTC  can  be 
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predicted  and  these  predictions  verified  [BIH,  1974). 

3. 4. 6.1  UH-UTC  Prediction 

BIH  tables  from  197'!  tv  i  /W  ’*»  V  >  e  analyzed  to  determine 
the  long-term  (>  6  month)  predictability  of  UT1-UTC.  Short-term 
predictions  have  been  made  by  Zhu  [  1 981  ]  and  Meyerhoff  [ 1 978 ]  in 
which  curve  fitting  techniques  were  used  to  fit  UT1-UTC  values  over 
one  year  (Meyerhoff)  and  three  years  (Zhu)  and  then  predict  these 
values  for  periods  of  5  to  40  days.  Each  study  used  a  series  of  the 
form 

UT1-UTC  -  a  ♦  bt  ♦  J  [Ci  sin  (  ^  )  ♦  C.  cos  (  )] 

(3.30) 

Zhu  set  n-2  and  Meyerhoff  determined  fits  for  n»1  to  15  with 
the  best  results  obtained  with  n»4  to  6.  Meyerhoff’ s  power 
spectral  analysis  indicated  that  n  must  be  at  least  4  to  fit  the 
major  frequencies.  Both  authors  then  predicted  UT1-UTC  values  for  a 
large  number  of  5-  to  40-day  intervals.  Meyerhoff  found  that  1o 

errors  of  2  to  7  ms  resulted  from  his  5-  to  20-day  predictions,  and 

Zhu  observed  prediction  errors  of  1.8  ms  for  5-day  predictions  to 

3.7  ms  for  40  days  but  that  the  errors  grew  rapidly  after  40  days. 

It  appears  that  the  determination  of  the  long-term  drift  over  three 
years  significantly  improved  Zhu's  results. 

Since  autonomous  operation  of  GPS  for  six  months  is  a 


goal,  the  long-term  prediction  of  UT1-UTC  is  a  requirement  of  the 
onboard  software.  To  determine  the  accuracy  of  this  prediction,  BIH 
data  from  1 97 ^  through  1980  were  fit  in  one-year  batches  by 
Eq.  (3.30)  with  n»4.  The  fit  resi;uals  shown  for  1974  and  1979 
(Figures  3.6  and  3.7)  are  typical,  and  it  is  seen  that  the  fits 
exhibited  maximum  errors  of  approximately  4  ms  and  RMS  values  of  2 
ms.  When  the  equations  were  used  to  predict  UT1 — UTC  for  one  year 
following  the  fit,  several  different  types  of  behavior  were 
observed.  The  best  prediction  was  for  1975  (Figure  3*8),  with 
maximum  errors  of  -10  and  +17  ms  and  an  apparent  long-term  periodic 
behavior.  1975  data  predicted  to  1976  (Figure  3.9),  however,  showed 
a  negative  slope  secular  trend  with  the  maximum  error  reaching  -78 
ms,  and  the  1979  fit-1980  prediction  (Figure  3.10)  had  a  positive 
slope  with  a  maximum  error  of  115  ms.  In  all  of  the  one-year 
predictions,  six-month  performance  was  better  than  70  ms. 


When 

data 

over  three 

years  were  fit, 

Eq.  (3.30) 

was 

augmented  by 

the 

addition  of 

two-year 

periodic 

terms. 

After 

correcting  for 

leap 

seconds ,  a 

sliding 

three-year 

fit  for 

data 

between  1974  and  1979  was  used  to  predict  UT1-UTC  values  for  the 
following  one  year.  The  results  for  one-year  predictions  at  six- 
month  intervals  are  shown  in  Figures  3.11  through  3«15,  where 
maximum  six-month  prediction  error  was  -58  ms  and  the  RMS  error  for 
all  predictions  was  36.7  ms.  Note  that  six-month  predictions  were 
much  better  than  those  for  one  year.  The  one-year  RMS  error  was  86 
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ms  and  several  examples  of  large  secular  growth  were  observed. 

From  these  data,  it  appears  that  errors  of  up  to  70  ms 
arise  when  one-year  BIH  data  is  used  to  predict  for  six  months  but 
that  three-year  data  spans  allow  predictions  with  1o  errors  of  less 
than  *40  ms  to  be  made.  Since  50  ms  appears  a  reasonable  limit  to 
six-month  UT1-UTC  predictions,  the  effect  of  this  error  on  GPS 
navigation  will  be  assessed. 

A  UT1-UTC  error  of  50  ms  corresponds  to  an  angular  error 
between  the  ECR  and  ECI  coordinate  frames  of  3.65  urad  *  .75  arcsee. 
At  GPS  altitude,  this  represents  an  apparent  ephemeris  error  at  the 
equator  of  96.5  m.  For  GPS  navigation,  this  is  a  large  error 
compared  to  the  desired  accuracy  but  may  not  be  as  severe  as  other 
error  sources  over  six-month  operation.  If  degraded  performance  is 
allowable  for  extended  autonomous  operation,  this  error  may  be 
acceptable.  In  event  this  error  must  be  reduced  better 
understanding  of  earth  rotation,  periodic  onboard  model  update  or 
active  GPS  tracking  of  earth-based  transmitters  would  be  required. 

3-^.6. 2  Polar  Motion  Prediction 

Both  Meyerhoff  [ 1 978]  and  2hu  [ 1981 ]  report  that  modeling 
of  the  moving  coordinates  of  the  geographic  north  pole  to  much 
better  levels  than  the  UT1-UTC  error.  Meyerhoff 's  20-day  pole 
positions  showed  x-only  coordinate  errors  of  approximately  .02  and 
.01  arcsec,  respectively,  while  Zhu  reports  60-day  average  RSS 
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errors  of  .02  arcsec.  Zhu  suggests  that  this  .02  arcsec  error 
remains  fairly  constant  for  up  to  two  years  of  prediction,  thus 
polar  motion  is  about  three  orders  of  magnitude  less  severe  than 
earth  rotation  rate  as  an  error  source. 

3. **.6. 3  GPS  Coordinate  Systems 

A  GPS  user  receives  satellite  position  information  in  the 
ECR  coordinate  system  [Van  Dierendonck,  et  al.,  1 980 ] ,  while  the  MCS 
software  determines  the  satellite  state  in  an  ECI  coordinate  system 
defined  by  mean  equator  and  equinox  of  1  January  1950  [General 
Dynamics,  1978]  and  uses  the  matrix  product  ABCD  to  transform  from 
ECI  to  ECR  coordinates.  D  is  a  rotation  matrix  containing  luni- 
solar  and  planetary  precession  terms  necessary  to  transform  from  the 
mean  equator  and  equinox  of  1  January  1950  to  the  mean  equator  and 
equinox  of  date.  C  is  a  matrix  containing  nutation  terms  necessary 
to  transform  from  the  mean  equator  and  equinox  of  date  to  the  true 
equator  and  equinox  of  date.  The  B  matrix  converts  from  the  true  of 
date  inertial  system  to  an  earth-fixed  system  by  rotating  through 
the  Greenwich  Hour  Angle  plus  UT1-UTC  and  nutation  (equation  of  the 
eqinoxes)  terms  and,  finally,  the  A  matrix  contains  polar  motion 
rotation  terms. 

The  onboard  estimation  algorithm  has  a  choice  of  these  two 
coordinate  systems  in  which  to  perform  its  ephemeris  calculations. 
If  the  ECR  frame  is  chosen,  conversions  of  output  are  unnecessary, 
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but  the  equations  of  motion  are  complicated  by 

“aj-Zw^v^ojxyxp 

where  <o  is  the  angular  velocity  of  the  earth.  The  errors  in  ui 
are  not  easily  modeled,  as  the  slope  of  the  UT1-UTC  curves  gives  the 
earth  rotation  error;  thus  errors  in  the  spacecraft  acceleration  are 
difficult  to  predict  in  the  ECR  frame.  In  addition,  the 
reconstruction  of  past  performance  of  the  navigation  filter  is 
complicated  by  the  time-varying  nature  of  6  w.  For  these  reasons, 
it  is  recommended  that  the  onboard  algorithm  use  an  inertial  (mean 
of  1950.0)  frame  for  the  ephemeris  calculations. 


CHAPTER  H 


GPS  ERROR  ANALYSIS 

*J.  1  Introduction 

A  relative  autonomous  navigation  application  for  PS 
spacecraft  will  include  an  observation  filter  and  an  onboard  model 
for  ephemeris  prediction.  This  model  of  the  satellite  dynamics  is 
used  for  propagating  the  spacecraft  state  and  error  covariance 
between  filter  updates.  The  propagation  interval  varies  from  the 
1.5  second  doppler  averaging  time  to  a  possible  one-  or  two-hour 
delay  between  updates. 

The  method  chosen  to  determine  the  requirements  for  the 
onboard  model  is  a  mixture  of  analytical  and  numerical  techniques. 
Satellit e-to-satellite  range  and  doppler  observations  from 
satellites  moving  in  similar  orbit  planes  cannot  provide  information 
concerning  a  common  secular  motion  of  all  satellite  planes,  but 
should  be  able  to  accurately  observe  differential  planar  and  in¬ 
plane  motions.  In  this  investigation,  secular  planar  motion  caused 
by  the  various  perturbing  forces  is  determined  analytically  through 
examination  of  equations  of  motion  of  the  orbit  elements  and 
numerically  through  a  consider  covariance  analysis,  while  the 
effects  of  errors  in  model  parameters  upon  planar  and  in-plane 
motion  is  determined  through  consider  analysis.  Figure  *1.1  depicts 
the  process  by  which  the  model  is  selected  and  analyzed.  The  model 


99 


101 


i 

l 


* 


I? 

n? 


9- 

k~ 


& 


4 


is  then  used  in  18  decentralized  estimation  filters  to  propagate  the 
state  vectors  of  the  18  satellites.  Observations  produced  by  a 
high-order  simulation  program  and  corrupted  by  random  measurement 
and  clock  noise  are  then  processed  by  the  estimation  algorithm.  The 
estimation  output  is  compared  to  a  truth  ephemeris  produced  by  the 
simulation  program  to  determine  the  accuracy  and  stability  of  the 
proposed  algorithm. 

Since  model  error  will  be  a  prime  consideration  in  the 
design  of  an  onboard  filter,  this  chapter  starts  with  a  discussion 
of  the  two  prevailing  techniques  used  to  handle  and  quantify  this 
error.  The  motion  of  GPS  orbit  planes  is  then  analyzed,  followed  by 
a  consider  analysis  of  the  expected  relative  navigation  accuracy  of 
the  system.  The  objective  of  this  chapter  is  to  define  a  model  to 
be  used  in  the  estimation  algorithm  proposed  in  Chapter  5,  along 
with  its  expected  performance. 

1) .2  Filter  Divergence 

As  derived  in  Chapter  2,  the  state  errc.  covariance 
matrix,  P  ,  when  propagated  and  updated  by  the  Kalman  filter 
equations  (Eq.  2. it  and  2.8]  gives  an  optimistic  (low)  value  for  the 
state  error  variance.  In  a  filter  where  the  state  is  being 
corrected  through  the  inclusion  of  observations,  this  low  value  of 
P  causes  a  low  Kalman  gain,  K  ,  which  results  in  new  data  being 
ignored.  If  the  true  spacecraft  state  is  given  by  the  n-vector 
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X(t)  and  an  a  priori  estimate  of  this  state  is  defined  by  X*(t), 
then  define  the  state  error  by 


x(t)  -  x(t)  -  x*(t)  KO 

The  nonlinear  orbit  determination  problem,  given  by  the  state 
differential  equation 

x(t)  -  F(x (t ) , t ]  (^-2) 


and  the  m-vector  measurement  equation  at  t  —  t ^ 

Yi  -  C(x(ti),ti)  ♦  e1  (4o) 

with  ei  -  n(o,R1),  Ri  >  0  , 

is  then  linearized  by  expressing  x(t)  as  a  Taylor  series  expanded 
about  the  a  priori  estimate,  X#(t)  and  Y^^  as  a  Taylor  series 
about  Y*  -  G(x*(ti ) ,ti ] .  The  linear  problem  is  then 

x(t)  -  A(t]x(t) 


yt  “  x(t1)  + 


where 


3G(x(ti),ti) 


axUi  j 


x* 


(tt) 


(^5) 


M 


An  unbiased,  minimum  variance  estimate  of  the  state  error 
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at  t  -  tk  is  then  given  by  Jazwinski  [ 1 970,  p.  198]  as 

*k  ‘  *k  +  KkK  '  Vk)  (4.7) 

where 

*k  “  ^k'VlVl  (4.8) 

and  ^(t4<*t({_1 )  is  defined  by  Eq.  (2.5). 

Note  that  in  Eq.  (4.?),  the  propagated  state  error 

estimate,  xk  is  updated  by  observation  residuals,  yk  ,  multiplied  by 
Kk  so  that  if  the  covariance  matrix,  "p  ,  drops  to  unrealistically 
low  values,  the  gain  Kk  -  pH^(hFH^  +  r)  1  becomes  small,  and  the 
updated  state  error  5?k  is  not  sensitive  to  the  observations.  This 
leads  to  the  well-known  symptom  of  filter  divergence  where  residuals 
become  large  and,  finally,  the  state  error  grows  to  values  much 

larger  than  the  error  bounds  expressed  by  the  variances  in  P  . 

The  cause  of  filter  divergence  is  the  assumption  that 
reality  is  correctly  modeled  by  the  two  sets  of  equations,  F(x(t),t) 
and  G(x(tj)ftj)  when,  in  fact,  our  knowledge  of  the  real  forces 
acting  on  a  spacecraft  and  of  the  actual  physics  involved  in  a 
measurement  process  is  quite  limited.  Even  if  we  had  perfect 
knowledge  of  the  forces  and  measurements  involved,  the  solutions  to 
the  state  differential  and  measurement  equations  are  only 

approximated  on  a  digital  computer.  The  two  major  errors, 
approximation  of  the  equations  solutions  by  finite  series  or 
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numerical  integration  (truncation  error)  and  finite  computer  word 
length  (rcund-off)  would  continue  to  give  errors  for  x(t)  and 

yi  • 

Filter  divergence  is  not  observed  in  pre-mission 
covariance  analyses  because  observations  are  not  being  processed  to 
produce  state  errors.  A  Monte  Carlo  analysis  in  which  a  more 

sophisticated  model  generates  observations  for  use  in  a  lower  order 
filter  can  be  used  to  evaluate  this  phenomenon.  Covariance  analysis 

results  do,  however,  suffer  from  the  effects  of  model  error.  The 

state  error  variance  is  optimistically  low  and  desired  information 
about  the  effect  of  unmodeled  errors  is  not  produced. 

4.3  State  Noise  Covariance  Analysis 

To  combat  the  problem  of  filter  divergence  or  an 

unrealistic  covariance  matrix,  estimation  algorithms  can  employ 
state  noise  or  model  noise  compensation  [jazwlnski,  1970,  pp.  244- 
247].  Instead  of  the  "perfect"  model  assumed  in  Eq.  (4.4),  let  us 
assume  that  our  state  differential  equations  are  corrupted  by  a 
Brownian  motion  process  such  that 

x(t)  «  A(t)x(t)  +  B(t)d8(t)  (4.9) 

-  Hix(ti)  +  el 

where  B(t)  is  a  non-random  matrix  that  maps  the  r-vector  d@(t) 
into  the  n  space  of  x(t)  and 


Ji  iVJf  J  i‘  1  11  r  -Twws*»»*#r/  ~*i*  ‘j-'-i'j.  ij,i»  *  *  ;• 


i; 


E[d8(t)]  -  0 

E[de(t)dgT(t)]  -  Q(t)dt  Q(t)  >  0 

E[dg(t )c16T(t) ]  -  Q6 (t**T )  6  ■  dirac  delta 

E[d8(t)e|]  -  0 


(*■•9) 

(4.<0) 

C**.n ) 

(**.12) 


The  state  error  is  then  propagated  from  tk  to  t  by 

,  ,  ,  ft 

x(t]  »  <Ht,tk)x(tk)  +  i)>(t , t)b( r)d8(t) . 


(*■•13) 


To  propagate  P  ,  first  look  at  the  propagation  of  P  without  model 


noise. 


F(0  -  ^(t.tk)p(tk)^T(t,tk) 


(**.1**) 


Differentiating, 


^(t]  -  ♦(tftk)P(tk)4t,tk) 

+  $(t »tk]p(tk)ij)  ( ) 
and  substituting  for  <j>(t,tk) 


(*■•15) 


”  A(t )<f>(t , tk ) 


(*■•16) 


we  arrive  at  the  matrix  Riccati  equation: 


F(t)  -  A(t)^(t,tk)P(tk)^T(t,tl<) 


+  ♦lt»tk)p(tk)*T(t.t|{)AT(t) 


-  A(t)F(t)  +  F(t)AT(t). 


Now,  for  the  case  of  model  noise,  if  8(t)  is  an  unbiased  estimator 
of  the  state  error  x(tj  and  t  >  tk  [Maybeck  1979,  pp.  16^1-167], 

X(t)  -  E[x(t)  |  yi  ,  i-i . «]  (4.18) 

and  the  covariance  is  defined  as 

P(t)  .  E[ (x(t)  -  «(ti)(x(t)  -  »(t))T] 

-  E[ax(t)ixT(t)]  (4.19] 


where 


Ax  -  x(t)  -  *(t)  «  4>(t,tk)x(tk) 


<J>( t , T )b( T )de C t)  -  ♦(t,tk_1  )s(tk_1 ) 


4>(t*tk)Ax(tk) 


kJ  + 


4>(t ,  T)B(T)d0(  t) 


then 


Mt)  -  E[  l$(t,tk)Ax(tk)  +  j  4>(t,T)B(T)d8(T)} 

tfck 

U(t,tk)Ax(tk)  +  |  4>( t ,  t)b(t)c1b(t)  }T] 

tk 

■  E[4>(t.tk)Ax(tk)AxT(tk)(fT(tk)] 


(4.20) 


‘  ’w* 

1  «- 

r  f 

\  ' 


+  E[*(t,tk)Ax(tk)  [  d0T( t)bT(t)^T^ ,  t  )  ] 


E[  ^  <|»(t , t)b( T)dB(t )  4>(t ,tk ) Ax (tk ) ] 

[|  f"  «^(t,T)B(T)d6(T)d6T(s)BT(s)^T(t,s)].  (H.2l) 


fck  fck 


Now,  assuming  that  Ax (tk)  and  dg(t)  for  tit. 


uncorrelated,  l.e., 


E[Ax(tk]d0T(t),  t  i  tk]  .  0 

and  noting  that  <J> ( t , t 3  and  g(t)  are  non-random,  and  that 


E[dfi(T)d8(s)]  -  Q6(t-s)  , 


[H.22) 


C^-23) 


F(t)  -  <f>(t.tk)p(tk)<(,T(t,tk) 


4l(t,T)B(T)Q(T)BT(T)QT(t,T)dT 


(4.24) 


Differentiating  by  Leibnitz  rule  and  using  Eq.  (4.16), 


F(t)  -  A(t)*(t,tk)p(tk)*T(t,tJ  +  <}>  ( t  *tlf)p(tLf)  d»T  Ct,tu)AT(t) 


k'  'vkJ 


+  $  ( t , t)b( t)q(t)bt( t ) ( t , t ) J “  l 

T  \ 

*  [A(t)«|>(t,T)B(T)Q(T)BT(T)d>T(t,T) 


♦(tft)B(T)Q(T)BT(T)#T(tiT)AT(t)]di 
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■F(t)  -  A(t)p(t)  +  F(t)AT(t)  +  B(t)Q(t)BT(t)  (^4.25) 

Eq.  I1!. 25)  shows  that,  since  Q(t)  >  0  ,  the  covariance  grows 
between  tk  and  tk+1  by  the  integral  of  B(t)Q(t)fiT(t)  ,  where 
the  matrix  Q  describes  the  error  variance  of  the  state 
differential  equation  model  error.  This  technique  of  augmenting  the 
P  matrix  by  a  Q  matrix  is  in  common  use  in  filter  algorithms.  The 
Space  Sextant  software  incorporates  the  Q  matrix  formulation 
[Martin-Marietta,  1 981 ]  and  it  is  used  in  the  GPS  navigation 
algorithm  proposed  in  Chapter  5. 

Sequential  Consider  Covariance  Formulation 

The  state  noise  compensation  method  described  in  Section 
14 - 3  prevents  filter  divergence  when  proper  values  for  state  noise 
are  used  in  the  Q  matrix.  However,  it  has  three  disadvantages 
when  used  for  covariance  analyses.  First,  the  values  describing 
state  differential  equation  noise  tend  to  be  arrived  at  in  an  ad  hoc 
fashion  since  the  actual  errors  are  unknown,  Just  as  the  actual 
force  acting  on  the  spacecraft  over  time  is  unknown.  Second,  and 
more  important  from  an  error  analysis  viewpoint,  is  the  difficulty 
in  isolating  and  evaluating  the  effect  of  specific  errors  in  the 
state  model  on  the  error  covariance.  Third,  it  is  impractical  to 
analyze  the  effect  of  estimating  a  parameter  that  is  causing  large 
errors,  because  the  entire  filter  must  then  be  redesigned. 

To  overcome  these  difficulties,  consider  analysis  has  come 
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into  widespread  use.  This  technique  allows  the  effect  of  errors  in 
model  parameters  upon  the  dynamic  3tate  to  be  considered.  Bierman 
[1977.  pp.  162-I7l]  describes  a  batch  processor  in  which  the  state 
vector  is  partitioned  into  a  set  of  variables  and  parameters  that 
are  to  be  estimated  and  a  set  of  parameters  whose  effect  is  only  to 
be  considered.  After  processing  a  fictitious  set  of  noisy 
observations,  his  algorithm  then  propagates  the  entire  covariance 
through  desired  time  intervals  to  determine  the  state  error  in  the 
future. 

In  the  GPS  scenario,  however,  observations  would  be 
processed  onboard  in  a  sequential  filter,  thus  the  sequential 
covariance  analysis  proposed  by  Maybeck  [l979,  pp.  325-336]  is  used. 
This  method  assumes  that  a  high-order  "truth  model"  is  available 
that  adequately  represents  the  real  world.  It  is  given  by  the 
linearized  continuous  differential  equation 


*t(t)  *  At(t)ct(t)  ♦  Gt(t)wt(t) 

h.28) 

• 

where  5t(t)  is  an  nt  vector  and  «t(t) 

mean  white  Gaussian  noise  sequence,  with 

is  an  St  vector  zero- 

E[et(t0)]  -  0 

(it. 29) 

B[«t{t0)eJ(t  )]  -  rt 

A 

(H.30) 

Discrete  measurements  available  from  this  model  are  obtained  by 


f 


* 


t 
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yt(tj)  ‘  +  etM  Ml) 

with 

E[et(ti)]  “  0  (4.32) 

^[^1  (^i  )et (tj  ) }  “  Rj.(t^) 

0  ti  -  tj  ,  (A. 33} 

For  consider  analysis,  the  state  vector  contains  both  the 

spacecraft  dynamic  state  and  constant  model  parameters.  The  dynamic 
state,  which  is  of  interest  at  each  time  t  ,  is  separated  by  £t(t) 
by 

xt(t)  -  Ct«t(t).  (4.34) 

After  arranging  in  a  manner  consistent  with  the  desired 

outcome,  the  matrix  can  be  partitioned  into  a  p  x  p  identity 

matrix  and  a  p  x  nt~p  null  matrix. 

The  "truth  model"  is  then  used  as  a  basis  with  which  to 
compare  several  reduced  order  candidate  linear  "filter  models,"  each 
described  by  the  n-vector  differential  equation 


l(t)  -  F(t)Cft) 


(4.35) 


with 


E[C(to)]  -  0 


(4.36) 


Ill 


-  rc  (1.37) 

and  updated  by  measurements  modeled  by 

?(*!)  -  H(ti)8Cti)  +  e(ti)  (K.38] 

with 

E[e(.tj)]  -  0  (4.39) 

E[e(ti)cT(tJ )]  -  R(tx)  |  tA  -  tj 

0  ]  *  tj  (4.  i40) 

The  Kalman  filter  for  these  filter  models  is  the  same  as  described 
in  Section  J). 2 .  Again,  only  a  subset  of  this  state  vector  contains 
dynamic  terms,  thus 

x(t )  -  Cf-(t)  (4.41 ) 


where  C  is  a  p  x  p  identity  matrix  plus  a  p  x  n-p  null  matrix. 

To  compare  these  two  models,  form  the  augmented  state 

vector 


(4.42) 


with 


pa(t0)  -  E[ea(t0)eJ(t0)] 


0 


(4.43) 


This  augmented  state  is  propagated  by 


“  *a^fck,fck-1^a^fck-J 


(4.44) 


where  <f>a (t^ * tk_i )  is  an  (nt+n)x^it+n)  state  transition  matrix 


that  satisfies 


^(t.t^.j)  -  Aa(  t  )<t>a  (t ,  tk_|  J 


(4.  45) 


^a^k-l^k-O  “  1 


(4.46) 


Aa(t) 


At(t)  0 
0  A(t) 


(4.47) 


Since  the  state  vectors  Ct(t)  and  |(t)  contain  both 
dynamic  and  constant  terms  and  the  consider  analysis  produces  a  time 
history  of  state  errors  due  to  these  terms,  it  is  convenient  to 
represent  <J>a  by  its  components.  The  two  dynamic  state  transition 
matrices  are 


3xt(tk) 


't^k-V 


(4.46) 


3x(tJ 


(4.49) 
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while  the  two  state  parameter  transition  matrices  are 


'Mtk»tk-1 ) 


3xt^k^ 


(1.50) 


3x(tk) 

3c,«v, ) 


(4.51 ) 


where  Cfc  and  C*  select  the  parameters  contained  in  5  and  £. 


The  augmented  state  covariance  is  moved  forward  in  time  vii 

^a^k^  “  ♦a^tk»tk-1^pa^tk-1^a^tk»tk-l^ 

Measurement  processing  gives  the  state  measurement  update; 

^(tk)  -  Ma(tk)Ta(tk) 


(4.52) 


(4.53} 


where 


MaUk) 


I 

_K(tk}Ht 


(tk)  fl-K(tk)H(t.<)  ] 


(4.54) 


and  Hfc(tk)  contains  partials  of  the  measurement  with  respect  to 
the  entire  truth  model  state. 

The  covariance  is  updated  by 

P.ltk>  ■ 


(4.55) 
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where 


Note  that  the  covariance  of  interest  is  actually 
Pjt)  -  E[nt(t)nJ(t)] 

where 

n(t)  -  x(t)  -  xt(t) 

This  covariance  is  obtained  from  Pg(t)  by 

Pn(E)  -  capa(t)cp 

where 


(4.56) 


b.57) 


(4.58) 


(^59) 


Ca  „  [-Ct  (  C]  (4.60) 

If  x(t)  consists  of  7  and  7,  then  P  contains  the  sum  of  the 

n 

6x6  upper  left  submatrices  of  Pt  and  P  . 

This  approach  to  state  covariance  analysis  assumes  that 

the  filter  operates  independently  of  the  higher  order  truth  model 

and  that  the  errors  in  the  truth  model  given  by  P,.  are  realistic. 

*0 

If  this  is  true,  P  (t)  will  represent  the  actual  errors  exhibited 
by  the  filter  in  use.  As  in  the  state  noise  compensation  technique, 
hovever,  the  actual  error  of  a  dynamic  or  geometric  parameter  may 
not  be  known.  In  this  study,  normalized  partials  are  used  to 
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determine  the  sensitivity  of  the  state  to  a  unit  error  in  each  of 
the  parameters. 

*•  • 5  Application  to  GPS  Autonomous  Uavlgatlon 

COVSEQ  is  a  computer  program  devised  to  implement  the 
sequential  consider  covariance  analysis  method  for  a  GPS  satellite 
navigation  receiver.  Variables  and  parameters  estimated  or 
considered  in  the  program  are  given  in  Table  <4.1. 

The  program  makes  use  of  the  propagation  capabilities  of 
the  University  of  Texas  Orbit  Processor  (UTOPIA)  to  obtain  a  user 
satellite  file  containing  the  spacecraft  state  (F(t),  v"(t))  and 

partial  derivatives,  <f>a ( t , tQ )  for  the  several  parameters  available 
in  UTOPIA.  Since  UTOPIA  does  not  produce  partials  for  the  state 
with  respect  to  lunar  and  solar  mass  and  position  errors,  these 
partials,  along  with  all  geometric  measurement  partials,  are 

computed  in  COVSEQ.  UTOPIA  produces  <f>( t , tQ )  and  ip(t , t0 ) ,  while 
COVSEQ  requires  and  <p(tj+1  ,tj ).  This  converrsion  is 

described  in  Appendix  A.  The  output  of  COVSEQ  is  a  time  history  of 
the  error  covariance,  P^  ,  after  a  set  of  GPS  range  and/or 
integrated  doppler  observations  is  processed  at  each  desired  step. 
The  error  covariance  is  propagated  according  to  Eq.  I1*. 52)  and 
updated  at  each  observation  time  by  Eq.  (i|.55). 


.{ 


5.1  Pseudo-Range  Observation  Partlals 

Measurements  available  to  the  GPS  vehicle  include  pseudo¬ 
range  and  integrated  doppler  counts.  For  the  crosslink,  the 
pseudo-range  model  is  given  in  Eq.  ( 3 - 1 )  and  state  and  consider 
partials  are: 


HxW 


x-x„  y-y_  2-zc 


ooo 


i  3pi 
C  3At 

u 


i  3pj 
c  9Afu 


1  8pl  1  3pl 
c  3At  c  3Af 

Si  SiJ 

(K.61) 


where  p  ,  At..  ,  and  At„  are  defined  in  Section  3.2  and  the  other 

A  U  3^ 

geometric  partials  are  given  by 


1  3pi 

— •  ■■  ■'  *  •  t  “  t 

c  3Afu  o 


(Jl.62) 


«  sensitivity  of  observation  to  a  user  clock  frequency  bias 


9 

3y7 


h.63) 


sensitivity  of  observation  to  GPS  ephemeris  errors 
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expressed  in  a  GPS-centered  RTN  coordinate  system 


(4.64) 


-  sensitivity  of  observation  to  a  GPS  clock  phase  bias 


3  Pi 


3  At 


c  3Af  "  "  lo)  (4.65) 

Si 

»  sensitivity  of  observation  to  a  GPS  clock  frequency 
bias 

^•5.2  Doppler  Measurement  Model 

When  a  signal  of  frequency  f  is  transmitted  by  a  moving 

O 

source  from  t1  to  t2  and  received  by  a  moving  receiver  between 
t1  and  t2  ,  the  received  frequency,  fr  ,  differs  from  fg  due  to 
the  well-known  doppler  shift.  When  the  received  frequency  is 
compared  to  a  known  oscillator  output,  fu  ,  the  resulting 
uifference  is  then  integrated  to  produce  the  doppler  count: 

fT  2 

N  "  J  (fu  -  fr)dT 
T1 

C2  C2 

-J  V1  "  J  frdT  »  (*».66) 

T1  T1 

but  fr  undergoes  the  same  number  of  oscillations  during  t2  -  t1 
as  does  f^  between  t2  and  t^  ,  thus 
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fudr  - 


dx 


Now  assuming  that  f  and  f 
u  g 

such  that  they  are  constant 
that  x  and  t  are  related  by 


M7) 

are  the  result  of  stable  oscillators 
during  the  intervals  used  and  noting 


x  -  t  +  £  (4.68) 

c 

where  p  is  the  geometric  range  and  c  the  speed  of  light 
(c  “  3.0  x  10®  m/s),  the  doppler  count  is  then 

N  "  fU(T2  "  *l)  "  fgK'  "  Tl)  +  *T  ^p2  "  Pi) 

-  (fu  -  fg)(x2  -  xt)  (P2  -  Pl)  (4.69) 

Since  Eq.  (4.69)  was  obtained  for  perfect  frequency  standards,  it 
must  be  modified  to  include  errors  in  frequency,  even  though  they 
are  small  for  GPS  clocks.  Expressing  the  true  oscillator  output  as 
a  nominal  frequency,  f*  ,  plus  a  first-order  term,  Af  ,  Eq.  (4.69) 
becomes 


N  “  ^fu  fg^T2  "  Ti  )  +  (Afu  “  Afg^T2  "  Tl) 


f  +  Af  _ 
8  8 


(P2  '  Pi  )■ 


(4.70) 


Noting  that  for  GPS-to-GPS  communication  f*  -  f*  •  and  multiplying 


both  sides  by 
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(*.71 ) 


Now,  f*  for  L3  is  approximately  109  Hz,  while  the  frequency 
stability,  Af/f  is  on  the  order  of  10-12,  thus  -j£-  <<  1  .  Note 
also  that  represents  a  delta  range  measurement,  thus 


Ap  -  c(«fu  -  «fg)Cx2  -  ♦  p2  -  p, 

where  6f  »  • 


(11.72) 


4.5.3  Doppler  Partial  Derivatives 

The  partial  derivatives  of  Ap  can  be  computed  at  either 

T1  or  t2  .  Assuming  that  the  receiver  integrates  from  t-j  to  t2 
and  then  estimates  the  spacecraft  state  at  t2  ,  the  observation- 
state  partial  derivatives  are: 

a.  Position  and  Velocity 

3Ap  _  3p(T2^  9p(T1 )  9X(Tl ) 

Pv  ^2 

-  — -  — L  — -  0  0  0 

P  P  P  Jt2 

“Pi  “  Pp  “Po 

-  _~-T  T00  °JX/W  (“-73) 

where  4> ( t1  , -r2 )  is  the  state  transition  matrix  from  x1  to  t2  . 


^  V»*SM  ^VVIi^rlSl 
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b.  User  Clock  Bias  and  Drift 


(4.74) 

u 

9Ap 

TOT  "  t2  '  T1 

(4.75) 

u 

GPS  Clock  Bias  and  Drift 

s 

(4.76) 

"  "(t2  "  Tl) 

(4.77) 

4.6  Model  Determination 

4.6.1  Geopotential  Coefficients 

The  major  geopotentiaal  source  of  secular  motion  o 
satellite  planes  is  the  well-known  J2  perturbation.  J2 
mathematically  models  the  equatorial  bulge  and  causes  a  regression 
of  the  nodes  given  by 


2 

n  cos  i 


where 

Re  -  equatorial  radius  of  the  earth 
a  «  orbit  semi-major  axis 
e  -  orbit  eccentricity 


(4.78) 


n  »  orbit  mean  motion 
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i  »  orbit  Inclination 


The  effect  of  neglecting  J,  is  sheen  in  Figure  «.2,  where  it  Is 

seen  that  the  error  normal  to  the  original  plane  of  873.6  km  builds 
up  over  60  days. 


The  other  main  geopotential  perturbations  are  due  to  terms 

WhiCh  ^  reS°nant  with  the  12^r  orbit  period.  These  terms  (i/2 
and  1A)  do  not  cause  large  short-term  (less  than  one  day)  out-of- 

Piane  motions  but  can  lead  to  along-track  errors  large  enough  to 
affect  user  navigation  accuracy.  Figures  4.3  through  4.8  show  the 
error  growth  due  to  neglecting  these  terms,  while  Figure  *1.9  shows 
orst  case  out  of-plane  error  due  to  neglecting  c.S  3/2.  The 
in-plane  error  plots  also  include  the  error  due  to  propagating 
initial  position  and  velocity  errors.  This  error  reaches  38. *159  km. 
and  when  it  is  removed  from  the  analysis,  the  remaining  in-plane  and 
cross-track  errors  are  as  shown  in  Table  *1.2. 


TABLE  *1.2 


CONTRIBUTIONS  DUE  TO  NEGLECTING 
GEOPOTENTIAL  TERMS 


Neglected 

Term 


2 

C7S  2/2 

C.S  3/2 
C.S  5/2 
C.S  7/2 
C.S  *1/11 
C.S  6/4 


Maximum  Error 
After  60  Days 
liSOlJ 


Contribution  of  Out-of-Plane 
Geopotentiai  Motion 

©  (m) 


206*1.3*10 

*1*1.669 
167.363 
39.038 
38.1(595 
53.32*19 
38. *1637 


1925.871 

6.21*1 

128.897 


0.569 
1  x  10~& 


21.086 

0.005 


873608.0 

17.6 

355.0 

0.0 

0.0 

0.0 

0.0 
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Figure  4.7 


Predicted  Position  Error  with  and  S 
Neglected. 


i 

J 


o 


Figure  4.9.  Predicted  Normal  Position  Error 
C^2  and  Neglected. 
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These  results  show  that  a  model  used  to  predict  over  time  intervals 
longer  than  one  hour  should  include  J2  and  C,S  3/2  and  4/4  terms, 
while  predictions  of  approximately  10  m  accuracy  over  one  hour  can 
be  achieved  with  J2  and  C,S  3/2.  Note  that  the  out-of-plane 
motion  due  to  the  J2~induced  regression  of  the  nodes  is  the  major 
secular  planar  effect  and  that  C,S  3/2  cause  a  smaller  but 
significant  motion  of  the  orbit  plane. 


4.6.2  Solar  Radiation  Pressure 

The  effect  of  neglecting  solar  radiation  pressure  is  shown 
in  Figure  4.10,  with  the  maximum  error  being  2606.9  m.  For  this 
case,  the  effects  are  almost  entirely  in-plane,  but  periodic  planar 
movement  can  be  determined  from  equations  as  given  by  Geyling  and 
Westerman  [l971,  p.  123]  as 

i  -  Fp  3in  a  3ln  io^sin  62  "  sin  ei)  („i79) 

n2  m  r0 


and 


ft  : 


Fp  sin  a  (cos  63  -  cos  61 ) 


n2  m  r 


(4.80) 


where  the  sun  is  assumed  to  be  in  the  equatorial  plane  at  an  angle 
a  from  the  satellite  ascending  node,  and 


F  -  4.74  x  10-6  (l  +  n)  A  newtons 


(4.81  ) 


POSITION  ERP.3R  SIGHfl  (METERS) 


for  a  satellite  of  mass  m  ,  cross-sectional  area  A  ,  reflectivity 

n  and  radius  rQ  that  enters  eclipse  at  8^  and  exists  at  82  . 

This  motion  only  occurs  when  vehicles  are  eclipsed  by  the 
earth,  so  the  onboard  model  must  include  eclipse  computations.  The 
fairly  large  in-plane  error  caused  by  neglecting  solar  radiation 
pressure  is  more  significant  than  a  similar  sized  geopotential 
effect,  because  solar  parameters  are  highly  variable,  whereas 
geopotential  terms  may  be  regarded  as  constant.  It  is  this 
variability  that  requires  the  solar  reflectivity  to  be  estimated  by 
the  onboard  filter. 


4.6.3  N-Body  Effects 

The  acceleration  of  a  satellite  near  the  earth  due  to  the 
presence  of  a  third  body  is  given  by  Geyling  and  Westerman  [l971,  p. 
113]  as 


(4.82) 


where  1^  and  refer  to  the  perturbing  body,  and  r  refers  to 
the  satellite  with  motion  referenced  to  the  earth.  When  transformed 
to  orbital  elements  and  averaged  over  one  satellite  orbit,  the 
results  show  a  periodic  effect  in  orbit  inclination  and  a  secular 
change  in  the  ascending  node  of 


3upr3 

- TT  COS 


i 


n 


rad/orbit 


(4.83) 
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For  GPS  orbits,  these  rates  due  to  the  sun  and  moon  are 

Aflgun  -  -5  x  10"6  rad/orbit  (4.84) 

Afl _  -  -1  x  10‘5  rad/orbit  (4.85) 

moon 

and  after  only  one  month  (60  orbits),  the  cross-track  position  error 
at  the  node  is  24  km. 

Note  that  Aflavg  ,  to  first  order,  depends  only  on  orbit 
inclination  and  radius;  thus,  each  GPS  satellite  will  experience  the 
same  effect,  i.e.,  the  whole  constellation  will  precess  and  the 
results  will  not  be  observable  from  satellite-to-satellite  tracking. 
Therefore,  lunar  and  solar  gravity  must  be  modeled  in  the  onboard 
propagation  software.  As  for  the  planets,  their  effects  are  at 
least  4  orders  of  magnitude  below  the  sun's  [Geyling  and  Westerman, 
1971,  p.  113];  thus,  6-month  errors  of  less  than  15  m  would  be  seen. 

4.7  Relative  Navigation  Accuracy 

To  predict  the  performance  of  a  model  when  used  in  a 
relative  navigation  filter,  the  consider  analysis  program  is  run  in 
a  filter  mode  in  which  observations  are  assumed  to  be  processed  from 
each  visible  satellite  every  60  minutes.  Model  errors  include  earth 
geopotential  terms,  as  given  in  Table  3.3>  and  solar  radiation 
pressure,  as  shown  in  Table  4.2.  This  table  also  shows  the 
observation,  satellite  ephemeris  and  clock  errors  simulated. 

When  errors  in  model  parameters,  other  satellite  positions 


and  clock  biases  are  considered  and  observations  processed,  the 
position  error  stays  below  7  m  for  60  days  (Figure  Ml).  However, 
as  shown  in  Figure  M2,  when  errors  due  to  a  clock  drift  of 
10  sec/sec  are  included,  the  error  grows  to  2.391  km.  Estimating 
both  bias  and  drift  along  with  the  vehicle  state  brings  the  maximum 
error  down  to  40.8  m  as  depicted  in  Figure  M3. 

**•8  Consider  Analysis  Summary 

The  model  suggested  for  the  GPS  autonomous  navigation 
algorithm  includes  lunar  and  solar  gravity  and  the  earth 
geopotential  through  and  C,S  3/2  and  requires  that  solar 
radiation  pressure  and  each  GPS  vehicle  clock  bias  and  drift  be 
estimated.  Including  these  terms  in  the  consider  analysis,  plus 
errors  described  in  Tables  3.2  and  M  for  relevant  terms,  gives  an 
overall  one-sigma  accuracy  of  40.8  meters  over  60  days. 


POSITION  ERROR  SIGMA  (METERS) 


TIME  (DATS) 

Figure  4.12.  Estimated  Position  Error  with  Clock  Bias 

and  Drift  Neglected  Using  Range  and  Doppler 
Observations. 


CHAPTER  5 


GPS  NAVIGATION  FILTER 

5. 1  Introduction 

The  onboard  navigation  filter  that  each  GPS  vehicle  uses 
to  determine  its  navigation  state  is  a  compromise  between  numerical 
accuracy,  size,  complexity  and  speed.  Since  current  microprocessors 
do  not  have  the  word  length  or  speed  of  mainframe  processors,  the 
filter  must  be  designed  to  fit  into  a  limited  storage  and  time 
environment.  As  in  the  previous  sections,  it  is  assumed  that  each 
visible  GPS  satellite  is  tracked  every  36  seconds  and  that  the 

i  pseudo-range  and  integrated  doppler  observations  are  smoothed  by  a 

- 

|  local  curve-fitting  procedure  to  produce  one  pair  of  observations 

I  per  satellite  each  60  minutes.  A  sequential  filter  then  processes 

I 

1  up  to  pairs  of  range  and  doppler  measurements  and  updates  the 

|  navigation  state  each  hour. 

An  alternate  approach  to  the  sequential  filter  is  to 
collect  smoothed  observations  for  several  hours  and  process  them  via 
a  batch  algorithm  at  the  end  of  the  observation  span.  This  requires 
the  covariance  matrix  and  state  to  be  propagated  over  much  longer 
intervals,  however,  and  would  require  more  accuracy  and 
sophistication  in  the  dynamic  model.  This  approach  is  worthy  of 
study,  but  this  analysis  assumes  a  form  of  the  extended  Kalman 
sequential  filter  (EKF )  to  minimize  propagation  intervals.  Other 
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studies  have  been  performed  on  satellite  navigation  algorithms, 
e.g.,  Tapley,  et.  al.  [l98l],  and  square  root  formulations  of  the 
EKF  have  proven  stable,  efficient  and  accurate  for  autonomous 
navigation  using  GPS  signals.  The  results  of  Tapley,  et  al.,  are 
directly  applicable  to  the  GPS  autonomous  navigation  problem  since 
they  studied  estimation  algorithms,  dynamic  models  and  numerical 
integrators  for  onboard  navigation  of  LANDSAT-D-type  satellites. 
The  main  differences  are  that  GPS  is  higher  (semi-major 
axis  «  26575  km  vs.  7087  km)  and  that  the  onboard  filter  is  actually 
a  local  filter  solving  a  part  of  a  global  problem.  The  local  filter 
aspect  is  addressed  in  Section  5.3. 

5.2  Filter  Model 

The  dynamic  model  selected  for  the  GPS  navigation  filter 
algorithm  differs  from  the  LANDSAT-D  model  in  that  the  higher 
altitude  puts  the  satellite  out  of  the  atmospheric  drag  regime  and 
into  the  area  where  solar  radiation  pressure  and  luni-solar  gravity 
terms  become  significant.  As  noted  in  Section  4.6. 3,  lunar  and 
solar  gravity  effects  are  unobservable  in  satellite-to-satellite 
tracking  if  the  satellites  are  in  the  same  orbits,  and  accurate 
models  for  these  effects  must  be  included  in  the  dynamic  model. 
Fortunately,  lunar  and  solar  ephemerides  can  be  predicted 
accurately,  and  their  masses  are  known  to  high  enough  accuracy  that 
their  gravitational  effects  can  be  modeled  instead  of  estimated  in 
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the  filter. 

Solar  radiation  pressure,  however,  is  not  as  well  modeled, 
and  the  satellite  coefficient  of  reflectivity,  n  ,  and  its  rate, 
n  ,  must  be  estimated  along  with  the  satellite  position,  velocity, 
clock  bias  and  drift,  as  discussed  in  Chapter  *J. 

The  state  vector  for  GPS  #i  is  then 


XH 
X5 

(5.1) 

x6 
x7 
x8 
x9 
X10 

5.2.1  Dynamic  Equations 

The  dynamic  model  used  to  propagate  the  state  by  the 
nonlinear  vector  differential  equation 

X(t)  ■  F(x(t),t) 


V 

CAt0 

CAf/f 

n 


is 
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where 

a'g  -  gravitational  acceleration  of  the  earth 
3g  «  acceleration  due  to  solar  radiation 
Fs  »  perturbing  acceleration  due  to  the  sun 
Fffl  -  perturbing  acceleration  due  to  the  moon 

and  the  stochastic  forcing  noises  have  the  following  statistics 

E[e(t]]£*o  (5.3) 

E[e(t)eT(t) =  q£6(t~r)  (5.4) 

where  &  denotes  the  3  x  1  vector  acceleration  process  noise, 
e"a  t  or  the  scalar  clock  frequency  drift  noise,  ec  . 

The  nonlinear  state  differential  equations  are  linearized 
about  a  nominal  state  vector  to  produce 

x(t)  -  A ( t )  x ( t )  (5.5) 


where 
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»<«>  •  i^feVil|x(t).x.(t) 


3v 

3v 

3v 

3v 

3v 

3r 

3v 

3cs 

3n 

3n 

3a 

3a 

3a 

3a 

3a 

3r 

3v 

3cs 

an 

an 

3c  r 

3cr 

3cr 

3cr 

3cr 
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3v 

3cs 

an 

3n 
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3ce 

3ce 
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3ce 

3r 

9v 

3 

3n 

an 

3n 

In 

an 

an 

an 

3r 

9v 

3cs 

an 

3n 

_  0 

0 

0 

0 

0  _ 

where  cr  is  the  clock  rate,  ce  is  the  clock  error  and  cs  is 
the  clock  state.  For  this  application, 


■43x3) 

1(3x3) 

0(3x2) 

0(3xl) 

0(3xl) 

a21(3x3) 

0(3x3) 

0(3x2) 

A2ij(  3x1 ) 

0(3xl) 

0(1x3) 

0(1x3) 

A33 ( 1 x2 ) 

0 

0 

0(3x3) 

0(3x3) 

0(3x2) 

0(3x1) 

The  entries  of  A  are  described  in  the  following  sections. 


5.2.2  Geopotential  Model 

As  discussed  in  Section 
with  the  exception  of  J2  and 
and  periodic  perturbations  that 


4.6.1,  earth  geopotential  terms, 

C,S  3/2  cause  along-track  secular 
are  small  during  the  one-hour 
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propagation  interval  and  are  observable  in  satellite-to-satellite 
tracking.  The  effects  of  J2  are  large  and  include  both  tangential 
and  normal  excursions  and  are  thus  unobservable  to  the  navigation 
system.  They  must  be  included  in  the  dynamic  model,  which  becomes 


ar 

e 


where 


-y  x,/r3  (l 


-y  x 2 /r  (l 


L"P  x3/r 


3 

3  (l 


a) 

а) 

б) 


r  -/  x2  ♦  x2  +  x2 
R2 

a  ■  1.5  J2  (l  -  5  x?/r2) 
r2  5 

R2 

B  -  1.5  j  _|  (3  -  5  x?/r2) 
r * 

J2  -  .001083 

For  the  dynamic  model  given  above, 


(5.7) 


(5.8) 

(5.9) 


(5.10) 

(5.11) 
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a13( 
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e 
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e 

a22 

e 

a23( 

13e 

*23e 

where 


-HUR!P 


mmmim 


U5 


‘l!  -  (l  -  3x2/r2)(l  ♦  a)  ♦  x2  (l  -  10x2/r2)  (5-13) 

e  rj  r'  0 


3J2Re 


12  “  ( 3ux  1  x2/r5 ) ( 1  +  a)  +  — ~~  x1x2(l  -  lOx^/r2)  (5.1^) 


3J2Re 


13  ■«  (3ux1x3/i  5)(l  +  a)  + — —■  x1  x3(6  -  10x2/r2)  (5.15) 


3j5r: 


a22  -  (1  *  3x|/r2)(l  ♦  a)  ♦  — ~  x|(l  -  10x2/r2)  (5.16) 


3J0R 


a23  *  (3ux2X3/r5)(l  +  a)  +  — x2x^(6  -  lOx^/r2)  (5.17) 


a33  “  -f  t1  “  3x2/r2)(l  ♦  g)  ♦  ^y£  x2(8  -  10x2/r2)  (5.18) 


3.  ,_Z\ 

I  U  •  ^ 

'e  rJ  "  r’ 

The  acceleration  due  to  C,S  3/2  is  computed  by  rotating 

the  satellite  position  vector  to  the  ECR  coordinate,  computing 

satellite  latitude  (<j>)  and  longitude  (i),  and  taking  the  gradient  of 

che  3/2  geopential  term: 


,R 


u32  “  r  ^'pr^3p32(sin  <j> ) 


(C32C0S  2X  +  S^?sin  2A ] 


32- 


(5.19) 


where  C^2  and  $32  are  the  spherical  harmonic  coefficients  and 
P32(sin  <p)  is  the  Legendre  Associated  Function  of  degree  3  and 
order  2  for  the  argument  sin  $  .  can  be  computed  recursively 


1 A  6 


from  the  equation 


P«m(sin  *)  -  P£.2|m  ♦  (21-1 )  cos  $  PJMf|1,_1  (5.20) 

where  Pj^  .  o  if  m  >  l  .  Thus,  P^sin  <j>)  is  obtained  from 

P32(sin4>)  -0  +  5  cos  <p  P22  (5.21 ) 

p22(sin<f>)  -  3  cos  <p  P^  (5.2 2) 

where 

P1l(sini}>)  ■  cos  4>  .  (5.23) 

Thus, 

P32 (si n<t> )  -  15  cos^  <t>  (5.2H) 


Entries  of  A21  fo r  C,S  3/2  are  not  computed  since  the 
e 

perturbation  due  to  these  terms  are  so  much  smaller  than  for  J2 
and  would  thus  affect  the  orbit  over  relatively  long  prediction 
times . 


5.2.3  Solar  Radiation  Model 

In  this  analysis,  the  same  solar  radiation  model  is  used 
in  the  onboard  filter  as  is  employed  in  UTOPIA.  The  acceleration 
dur  *-o  radiation  pressure  is  given  by  McMillan  [  1 97 3 ]  as 


a 

s 


Ps  rs  ^ 


+ 


y 


(5.25) 
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where 

ps  “  solar  radiation  pressure  on  a  black  body 
rs  *  sun’s  distance  from  earth 
h  -  satellite  surface  reflectivity 

A  satellite  effective  area 
m  -  satellite  mass 

^vs  “  sun-satellite  radius  vector  «  r-7 

s 

y  -  shadow  indicator 

For  a  spherical  earth  and  cylindrical  shadow, 

Y  -  0  if  F  ‘  ~  <  o  and  (r2  -  r2)1/2  <  r 

Y  «  1  otherwise 

The  partial  derivative  matrix  is 


a1 1  al2  a 


-J 


S  '-s  Js 
'2,s  "  |  X  a223  a23s 


a13s  a23s  a33s 


(5.26) 

(5.27) 


(5.28) 


where 


a11  -  P3rf(l^n)  J±-  1-3. 


(x-xs)‘ 


-3p  r^..np  U-3H^a) 

SV  Um  c: 

rJ 

vs 


(5.29) 

(5.30) 


2H  is  then  a  vector  given  by 


Ar 


psrs 


vs 


mr 


vs 


(5.35) 


5.2. 4  Lunar-Solar  Gravity 

The  perturbations  due  to  lunar  and  solar  gravity  are  given 


by 


where 


^  - 


(5.36) 


vector  from  earth  to  sun 
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GM„ 


and 


where 

-  vector  from  earth  to  moon 

~K  -  F  -  F 


(5.37) 


ym  »  CMm 

The  vectors  Fg  and  7*m  are  computed  by  assuming  Keplerian  orbits 

for  the  earth  and  moon  and  solving  Kepler's  equation  for  each  time 

point.  The  partial  derivatives  of  V  and  7  can  Then  be 

s  m 

obtained  by  a  straightforward  differentiation  of  Eqs.  (5.35)  and 
(5.36). 


5.3  Decentralized  Filtering 

The  problem  of  estimating  the  state  vector  of  each  of  18 
GPS  satellites  from  satellite-to-satellite  pseudo-range  and  doppler 
information  is  a  global  estimation  problem,  since  the  measurement 
errors  are  a  function  of  other  satellite  position  errors  as  well  as 
ranging  system  errors.  To  solve  this  problem  on  a  global  scale 
would  require  a  "supervisor"  system  with  knowledge  of  all  spacecraft 
and  measurement  errors.  If,  at  the  other  extreme,  the  situation  is 


4  k 
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handled  by  18  decentralized  processors,  each  solving  the  local 
estimation  problem  with  no  cognizance  of  the  other  spacecraft 
errors,  an  unstable  solution  results.  Any  GPS  position  error 
outside  the  range  expected  by  the  a  priori  measurement  covariance, 
R  ,  would  quickly  lead  to  filter  divergence  in  the  other  processors. 

This  situation  is  similar  to  that  faced  by  the  Joint 
Tactical  Information  Dissemination  System  (JTIDS)  in  its  aircraft 
relative  navigation  mode  as  described  by  Kerr  and  Chin  [ 1 980 ] . 
Since  the  two  problems  are  alike  in  that  each  vehicle's  state  is 
independent  of  the  others  but  measurements  involve  information 
exchange  among  the  members,  it  is  possible  to  cast  GPS  autonomous 
navigation  as  a  relative  navigation  problem  and  apply  known 
techniques  of  decentralized  filtering. 

Applying  the  methods  described  by  Kerr  and  Chin  to  the  GPS 
problem,  the  system  state  error  is  described  by  a  linear 


differential  equation  of  dimension  18  n^  x  i 

dimension  of  each  spacecraft  state 

where 

n^  is  the 

A(t)  -  F ( t )  : 

x(t)  ♦ 

w(tj. 

(5.38) 

This  system 

can 

be  expressed  as 

the 

collection 

jSi  .  i-1,  2,.. 

.,18  J 

of  separate  dynamic  systems: 

Si  :  XA(t)  - 

F±(t ) 

x±(t)  +  wi(t) 

(5.39) 

having  interconnected  discrete  measurements  available  to  the  ith 
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satellite: 


where  Mi(t|<)  is  the  projection  operator  from  Rn  to  R  1  ,  i.e. 

xi(tk)  -  Mi(tk]  £(tk)  (5.H1) 

and  Li  )  selects  those  components  of  X.(tk)  that  appear  in 

8G(x(tk,tk)/3Xi(tk)  .  Thus 

yi(tk)  -  H ( t k )x± (tk )  +  Hi(tk)Li(tk)x(tk)  +  vi(tk)  (5.42) 

and  contains  no  component  of  x^(tk)  directly. 

It  is  assumed  that  w^t)  and  vjtj  are  zero-mean  white  noise 
processes  that  are  uncorrelated  with  Wjft),  Vj(t)  and  x.(o)  for 
j  *  i  . 

Looking  at  the  specific  example  of  GPS-GPS  pseudo-range, 
assume  that  vehicle  i  is  tracking  vehicles  1-3  at  t  -  tk.  The 
measurement  matrix,  1L  for  a  satellite  state  vector  of  dimension 
10,  containing  r,  v,  clock  bias,  clock  drift,  reflectivity  and  its 


rate  is 


152 


for 


V 

\ 

V 

\ 

x  -  X 

i  J1 

pl 

pl 

pl 

\  “ 

*1 

X2  " 

x2„ 

X  -  X 

i 

2 

i 

2 

i  J2 

P2 

P2 

P2 

*1.  " 

xl. 

X2 .  ' 

*2, 

x3  "  x3 

1 

3 

i 

3 

i  3 

p3 

P3 

P3 

0  0  0  1  0  0  0 


0  0  0  1  0  0  0 


0  0  0  1  0  0  0 


(5. *3) 


r  v  cAto  n  n 
1  o 


(5.  *1*0 


5.3.1  Decentralized  Filter  Algorithms 

As  discussed  by  Kerr  and  Chin,  the  Surely  Locally  Unbiased 
(SLU)  [Sanders,  et  al . ,  1973]  and  the  Sequentially  Partitioned 
Algorithm  (SPA)  devised  by  Shah  [ 1 97 1 ]  are  stable  solutions  to  the 
relative  navigation  problem.  While  the  SLU  filter  is  unbiased,  it 
requires  that  the  rank  of  H  be  less  than  the  row  dimension.  This 
is  not  satisfied  for  the  pseudo-range  observation  [and  for 
integrated  doppler),  thus  the  SLU  algorithm  is  not  applicable. 

According  to  Kerr  and  Chin,  the  SPA  filter  has  been  shown 
to  be  asymptotically  stable,  but  it  is  not  analytically  unbiased. 
They  state,  however,  that  this  potential  problem  can  be  handled 
adequately  and  recommend  the  SPA  as  the  algorithm  for  JTIDS 
navigation  processors.  Biased  solutions  are,  however,  a  potential 
problem  area  in  GPS  navigation  and  must  be  investigated. 
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5.3.2  Derivation  of  the  Sequentially  Partitioned  Algorithm 

Following  Kerr  and  Chin  [ 1 981 ]  with  appropriate 
modifications  for  GPS,  define  the  local  system 


xi(k+1  ) 

-  4>i  (k+1  ,k)x1(k)  +  w^  (k ) 

(5.^5) 

with  observations 

y±(«<)  - 

H(k)x(k)  +  vi(k) 

18 

where 

Hi(k)xi(k)  +  Hjk)  l  Lij(k)xj(k)  +  vi(k). 

(5.^6 ) 

■  p  '  nj  x  1  vector  for  p  satellites  with  n^ 

each 

states 

«  n^  x  1  state  vector  for  satellite  i 

Lij 

-  m  x  nj  submatrix  relating  m  observations 

satellite  j  to  satellite  i 

3G(x(tkJ,tk) 

3x(tk) 

from 

H(k) 

(5.'«7) 

3G(x(tk],tk) 

‘“iW 

(5.H8) 

Hi(k) 

M*) 

3G(x(tk),tk) 

(5.«9) 

8U(tk)'xi(tk)l 

Now,  define  the  state  estimation  errors 


15*. 


ei(l<|k)  =  x.(k)  -  ^(k)  (5.50) 

and 


(k | k— 1 )  =  x.(k)  -  x.(k) 

(5.51) 

where 

Ck )  »  e[x± (k) |y(k) ] 

(5.52) 

and 

*iw  ■  nivW 

(5.53) 

Then 

-  a  18 

yi(k)  -  H.(k)x  (k)  *  H  (k)  I  L  (k)7  (k)  ♦  v*(k) 

j=1  J  1 

(5.5*l) 

where 

* .  .  1 8 
vjk)  =  v.(k)  +  H.(k)  l  L  (k)e.(k|k-l) 

j  =  1  J  J 

(5.55) 

Now,  assuming 


1.  Xj(k)  is  known  for  j*i,  j=l,  2,  ...,  18 

2.  e j C k j k ) ,  ej (kjk-1 )  are  Gaussian  and  white 
J-1,2,...,18,  j*i  , 
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In  the  extended  filter  where  the  state  correction  is  added 
to  the  state  after  each  measurement  update, 

xjk+l)  =  X.(k+l)  +  5ci(k+l)  (5.6l) 

where 

Xjk+l)  -  K.(k+l)yi(k+l)  .  (5.62) 

Comparing  the  SPA  to  consider  covariance  analysis,  it  is 
seen  that,  for  no  interconnection  in  the  dynamic  model,  the  SPA  is  a 
consider  filter  in  which  errors  in  GPS  vehicles'  positions  and 
clocks  are  considered. 


5.3.3  Sequential  Processing  of  Pseudo-Range  Observations 

For  the  case  of  GPS  i  obtaining  pseudo-range 
observations  from  one  satellite  (j)  at  a  time  in  a  sequential 
receiver,  assuming 


*fi  * 
r  v  At.  — 


Hi(k) 


“(x1  -xi  )  -(: 


tx  -x  )  -(x  -X  ) 

j  l  j  i  3  3  1 


000-1  000... 


. .  .0  0. . . 


(X1.“XJ  (x2.-x2)  (x3."X3.] 

j _ 1  J  1  J  1 


000-1000. 


(5.63) 
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with  zeroes  for  all  partials  with  respect  to  X&  ,  i  *  i  or  j  . 
This  vector  is  then  separated  into 


H1(k) 


{1  "xi  )  “( 

J  i 


x  -x  )  -(x  -X,  . 

J  j  3i 


0  0  0  -10  0  0 


=  [  h1  h2  h3  0  0  0  h?  0  0  0  ] 
and 

A 

H± (k)  -  [  -h1  -h2  -h3  0  0  0  -h  0  0  0  ] 
The  L  matrix  for  this  situation  is 
Lu  *  [  0  ]  for  i  *  j 
and 


(5.64) 


(5.65) 


(5.66) 


(5.67) 


Now,  define 


« ■  «,(>■)  Ltt  pt  jit 

-!h  ,  h2  n3  0  0  0  h7  0  0  0]  LjjPjL^lh,  h2  h3  0  0  0  ;7  0  0  0]T 
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When  processing  integrated  doppler  measurements  in  a 
sequentially  partitioned  algorithm,  the  measurement  matrix  is 


4>(  T1  ,t2) 
(5.70) 


where  t  is  the  receiver  time  tag  and  x^  =  tk  .  For  the  short 
time  interval,  t2't1  =  1*5  sec  ,  thus 

♦  :  I  (5.71 ) 

and 

(k)  :  [(h1 (x2)-h1 (t1 ))  (h2(x2)-h2(T1 ) 

(h3(x2)-h3(T1  ) )  0  0  0  0  h8  0  0  0. . . (h1 (t1 )-h1 (x2) ) 
(h2(t1)-h2(T2)(h3(x1)-h3(T2))  0  0  0  0  -hg  0  0...](5.72) 
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where 


h8  *  t2  "  T1 


(5.73) 


As  in  the  case  with  pseudo-range  observations,  this  vector 
(for  one  observation)  can  be  separated  into 

Hjk)  =  [(h1(i0)-h1(T1))  (h2(T2)-h2(i1  )) 


(h3(T2)-h3(t1 ))  0  0  0  0  hg  0  0  ] 


(5. 7*0 


H.(k)  »  (-(h1(‘t2)-h1(t1))  -  (h2(t2)-  h  2(t1)) 


-(h3(T2)-h3(t1))  0  0  0  0  -hg  0  0  ] 


(5.75) 


Lu  «  [  0  ]  1*  j 


(5.76) 


L.  .  =  1  . 

ij  1 


(5.77) 


Now,  define 


B  =  H.(k) 


^  Li£P£Li£ 
£*i 


(5.78) 


ft 


...IJ.  L.  IU1JU.I.I,!IUH. 


•* 
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and 


Ah.  =  h.(t2)  -  h.ltj) 


(5.79) 


Then 


6  •  A^P,,  .  ^  *  h|p88 

J  J  J  j 


+  2(Ah1Ah2P12  +  Ah1Ah2P 

J  J  j 


+  Ah2Ah3P23j 


+  Ahlh8Pl8j  +  Ah2h8P28j  +  Ah3h8P38.) 


(5.80) 


and 


R* (k+1 )  =  R.(k+l)  ♦  e  . 


(5.81 ) 


Note  that  the  transmission  of  data  from  satellite  j  now 


requires 


p  p  p 

11  12  13 

P17 

P18 

P22  P23 

P27 

P28 

?33 

P37 

CO 

,  00 

F77 

P88 

(5.82) 


5. 4  Numerical  Results 


I 


The  Sequentially  Partitioned  Algorithm  was  included  in  a 
computer  program  (GPSNAV)  that  solves  for  corrections  to  the  state 


itriffiffli 
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vector  defined  by  Eq.  (5.1 )  for  each  of  the  18  GPS  spacecraft.  The 
original  program  employed  a  U-D  filter  as  described  in  Appendix  B, 
however,  numerical  instability  in  the  propagation  of  the  P  matrix 
on  the  Cyber  computer  led  to  the  decision  to  use  a  standard  EKF 
formulation  with  Joseph's  form  of  covariance  update  and  propagation 
via  the  state  transition  matrix. 

Each  spacecraft  receives  range  and  doppler  from  all  other 
visible  satellites  each  60  minutes.  It  also  receives  the  covarience 
information  defined  in  Eq.  (5.82)  at  the  time  of  transmission  from 
each  of  the  other  vehicles.  Range  and  doppler  Information  was 
generated  from  a  UTOPIA  model  incorporating  an  6x8  earth 
geopotential,  solar  and  lunar  gravity,  solar  radiation  pressure  and 
clock  bias  and  drift  errors  using  the  force  models  discussed  in 
Chapter  with  parameter  and  initial  state  errors  as  given  in 


Table  5.1 . 


Range  errors  of  1o  =  2.0  m  and  doppler  errors  of  1 o  = 


.001  m/s  were 
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added. 

The  model  incorporated  in  the  SPA  filter  is  a  reduced 
order  geopotential  consisting  of  and  C  and  S  3/2  plus  solar  and 
lunar  gravity  and  solar  radiation  pressure.  The  filter  estimates 
all  of  the  terms  in  Eq.  (5.1)  and  has  the  capability  to  include 
model  noise  (q  matrix)  terms  for  acceleration  and  clock  errors. 

The  program  was  run  in  a  perfect  clock  mode  in  which  clock 
errors  were  neither  added  to  the  data  nor  estimated  by  the  filter. 
The  results  for  satellite  //I  (Fig.  5.l)  show  that  the  corrections  to 
r  and  the  position  variance  stay  below  5  m,  while  the  observation 
residuals  are  on  the  order  of  10-15  m  and  that  the  solution  appears 
stable  for  the  four-day  run. 

When  clock  errors  are  included  in  the  data,  however,  the 

results  in  Figure  5.2  show  unstable  results  after  about  80  hours. 

This  run  includes  no  model  noise  for  clock  parameters.  When  a  clock 

-12 

drift  noise  term  of  10  m/s  is  included  in  the  noise  compensation 
matrix  (Fig.  5.3)»  the  solution  is  stable  for  approximately  17  days 
with  corrections  on  the  order  of  5  m  and  observation  residuals  of 
10  m  being  seen. 

When  a  four-day  filter-determined  ephemeris  is  compared 
with  the  UTOPIA  ephemeris,  it  is  seen  that  the  in-plane  components 
agree  to  within  40  m  (Figs.  5.4  and  5.5)  with  a  small  secular  trend 
(2  m/day)  apparent  in  the  along-track  direction.  Out-of-plane 
errors,  however,  grow  secularly  at  approximately  20  m/day.  This 
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Figure  5. A.  Radial  Position  Difference  Between  Estimator 
and  UTOPIA 
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difference  may  be  due  to  a  combination  of  integrator/model 
differences  between  UTOPIA  and  errors,  an  inability  of  the  filter  to 
determine  adequately  the  cross-track  component  of  the  vehicle's 
clock  error,  or  biasing  of  the  filter  algorithm.  When  clock  errors 
are  removed  from  the  data  and  the  filters  assume  perfect  clocks,  the 
error  growth  is  reduced  to  11  m/day  (Fig.  5.6),  thus  almost  50 
percent  of  the  error  is  due  to  the  clocks.  Since  neither  of  the 
observable  components  (radial  and  tangential)  show  significant 
biases,  this  very  limited  test  indicates  that  the  problem  discussed 
in  Section  5.3.1  may  not  be  significant. 

If  this  secular  trend  continues,  the  total  cross-track 
error  would  grow  to  1.2  km  after  60  days,  but  only  a  small 
percentage  of  the  error  would  appear  in  the  user-GPS  line  of  sight. 
However,  the  importance  of  this  error  growth  is  not  its  magnitude 
but  the  fact  that  any  onboard  orbit  determination  scheme  based  upon 
satellite-to-satellite  tracking  will  not  observe  similar  planar 
motion.  To  solve  this  problem,  either  the  model  must  be  tuned  to 
the  specific  application  through  much  more  extensive  analysis  and 
testing,  or  the  navigation  system  must  be  augmented  by  a  GPS 
transmitter  in  a  different  altitude  orbit  or  on  the  surface  of  the 
earth  or  moon.  A  transmitter  on  the  earth  would  solve  this  problem, 
plus  it  would  allow  the  navigation  algorithm  to  determine  the 
earth-satellite  orientation. 
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CHAPTER  6 


CONCLUSIONS  AND  RECOMMENDATIONS 


6.1  Conclusions 

The  conclusions  drawn  from  the  discussion  in  the  previous 
sections  can  be  summarized  as  follows. 

1.  When  compared  with  the  space  sei:tant,  current  design  solid- 
state  matrix  star  sensors  coupled  with  current  horizon  sensors 
yield  a  factor  of  16  worse  navigation  accuracy. 

2.  If  matrix  sensors  can  be  used  to  measure  the  earth  horizon  by 
star  refraction  to  a  precision  comparable  to  the  star  position 
measurement,  they  can  provide  better  navigation  accuracies 
than  the  space  sextant.  A  critical  factor  in  this  scheme  is 
the  fact  that,  for  a  constant  tangent  height  error,  apparent 
horizon  sensor  error  quickly  drops  to  levels  near  or  below  the 
original  actual  sensor  error  ai  altitude  increases.  A  key 
factor  limiting  refraction  determination  of  earth  horizon  is 
the  accuracy  to  which  atmospheric  density  can  be  predicted. 

3.  Since  matrix  sensors  can  be  derated  to  match  less  stringent 
program  requirements,  they  can  be  used  to  provide  angular 
information  with  a  precision  varying  from  >20  arcsec  to 
<1  arcsec  with  corresponding  navigation  accuracy. 


173 

The  space  sextan',  has  a  greater  ability  to  recover  from  loss 
of  attitude  control  than  do  the  fixeu  sensors. 

Atmospheric  refraction  uncertainties  contribute  approximately 
1500  m  to  the  error  in  the  determination  of  ray  tangency 
altitude. 

For  navigation  requirements  of  50  m  or  less,  the  GPS  receiver 
is  recommended.  It  is  possible  that  a  SHAD  and/or  a  matrix 
sensor  may  produce  real  accuracies  of  <50  m,  but  an  analysis 
of  all  relevant  error  sources  is  required  to  verify  this 
conjecture. 

An  onboard  navigation  system  for  GPS  based  upon  satellite-uj- 
satellite  data  transmission  is  feasible.  However,  this  system 
is  inherently  unable  to  determine  similar  motion  of  all 
orbital  planes  and  the  earth.  While  earth  angular  position 
can  be  predicted  to  about  50  ms  over  six  months,  this  error, 
in  addition  to  planar  motion,  would  cause  errors  normal  to  the 
satellite  plane  to  reach  approximately  1.3  km. 

The  Sequentially  Partitioned  Algorithm  provides  a  stable 
solution  for  each  satellite's  state  vector,  including  clock 
parameters,  as  long  as  model  noise  compensation  is  used. 
Preliminary  results  do  not  show  a  tendency  toward  biasing. 
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6.2  Recommendations 

These  recommendations  are  made  following  the  previous 

study. 


1.  Development  of  the  matrix  sensor  technology  should  be 
continued. 

2.  Investigation  of  the  feasibility  of  performing  improved 
horizon  determination  from  star  refraction  measurements  should 
parallel  sensor  development.  Current  and  future  refraction 
data  should  be  used  to  improve  atmospheric  density  modeling. 
This  is  critical  to  the  success  of  navigation  by  matrix 
sensors. 

3.  Further  refinements  should  be  made  in  the  analysis  of  the 
matrix  sensor  measurement  to  include  all  error  sources 
affecting  the  measurement  and  its  reduction  to  a  navigation 
measurement. 

4.  Further  studies  of  GPS  cross-link  navigation  should  be 
conducted  to  include: 

a.  A  more  detailed  evaluation  of  filter  performance  with 
respect  to  model  and  clock  parameters. 

b.  Evaluation  of  navigation  performance  when  augmented  by 
an  additional  transmitter  in  space,  on  earth  and  on  the 
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moon . 

c.  The  feasibility  of  including  navigation  information  in 
the  cross-link. 

d.  U-D  filter  performance  when  used  with  the  SPA  algorithm. 


APPENDIX  A 


CONVERSION  OF  BATCH  PARTIAL  DERIVATIVES  TO  SEQUENTIAL 

The  University  of  Texas  Orbit  Processor  (UTOPIA)  produces 
precise  estimates  of  position  and  velocity  for  spacecraft  using 
several  types  of  observations,  including  laser  range  and  altimeter 
height.  The  basic  batch  and  sequential  algorithms  are  described  by 
McMillan  [1973]  and  Wilson  [ 1976]. 

A  relatively  recent  addition  to  UTOPIA  is  the  capability 
of  producing  partials  of  the  spacecraft  state  at  observation  times 
(inertial  r(t^J  and  v(tj))  with  respect  to  the  state  and  several 
model  parameters  at  some  user-defined  epoch.  These  model  parameters 
include  ,  the  spherical  expansion  coefficients  (clm,  Slm), 
atmospheric  drag,  solar  radiation  pressure  and  unknown  spacecraft 
accelerations  in  the  RTN  coordinate  system.  The  purpose  of  the 
partials  is  to  allow  the  program  to  interface  with  the  JPL  consider 
covariance  analysis  program  COVAN  to  produce  a  consider  covariance 
matrix  at  epoch  (tQJ  after  processing  simulated  observations  at 
i»1,...,k.  UTOPIA  does  not  simulate  satellite-to-satellite 
range  or  doppler  observations  nor  does  it  produce  a  sequential  form 
of  the  $  and  \p  matrices  required  for  long-term  sequential 
analysis  of  the  GPS  autonomous  navigation  problem. 


It  is  possible,  however,  to  convert  the  batch  or  epoch 
forms  of  these  matrices  to  a  sequential  form  by  the  following 
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algorithm: 


A>1  Conversion  from  4>(tj+1  ,tQ)  to  ♦  (t^1,tj 

Given  the  sequences  *(t0,tj,  ♦(t1,t0). 

♦^2.t0)...*(tJ+1,t0)  and  *(t0.t0).  ♦(t,  ,t0).  *(t2,t0)...*(tJ+1#t0), 

we  desire  the  sequences  ♦(t0it0),  ♦(t1.t0)#  *(t2.t, ). . .*(tj+1 ,tj ) 

o*^o )  ’  'l'  C 1 2  *  t  -j  ) . . .  ij»(tj  +f ,  tj  ) .  Suppose  a  state 


is  composed  of  variables  X 

and  parameters  Z  ,  then 

ax(tj 

-1 J  "  , 

ax(t .) 

J  i 

(A.  1  ) 

axCtA) 

J  az(tj) 

5* 

(a. 2) 

Note  that 


3X(t0)3x(tj) 

so,  for  i  «  j+1,  the  sequential  form  of  <p  is 
-  ♦(tj+1,t0)  ♦(t0,tj) 

■  *(Vl’t0)  f'l'i'to)-  (A. 3) 

This  conversion  of  <j>  from  epoch  to  sequential  requires  the 
inversion  of  a  6  x  6  matrix  at  each  step.  This  can  be  avoided  when 
the  spacecraft  force  is  conservative,  and  thus  only  a  function  of 


position,  by  using  the  symplectic  property  of  the  A  matrix. 
Recalling  that 
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3F 

3r 


r3  r5  j”1! 

3Mx1X2  ,  ^2 
r5  +9x1 

3yxlx3  3P3 
r5  +3xl 


3yxix2  ,  9P1 

r5  9X2 

SW 2  9F2 

r3  r5  +  3x2 

3yx2x3  9P3 

5  +3x, 

r  2 


3yx^x3  3Pj 

3yx2x3  9P2 

5  3  x- 

r  3 

,  3Ux^  3P_ 

ZM  + _ 1  3 

3  5  SxT 

r  r  3 


(A. 7) 


and,  if  P  is  a  conservative  force,  then 


7  x  P 


(A. 8) 


where 


7 


k. 


For  Eq.  (A. 8)  to  hold, 


3P1  ^  3P2 
'5x'2  ”  &x1 

3Pt  3P3 
3x'3  “  3x^ 

3P2  3P3 

3x3  ”  3x2 


(A. 9) 


and  the  submatrix  iL  is  symmetric. 

3F 

The  symplectic  property  of  A(t)  then  states  that  if  A ( t ) 


can  be  expressed  as 
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11 

A1 2 

'21 

A22 

with 

A  ,T 

*11  *  “A22 

A  J 

A1 2  *  a12 

A21  "  AIl 

then 


(a.io) 


(a.h) 


4»(t  ,t0) 


$11  $12 
$21  $22 


and 


(:  12) 


T  T 
$22  ”$1 2 
,T  <fcT 

-$2i  $11 


(A. 13) 


As  discussed  in  Goldstein  [ 1 980 ] ,  the  state  transition  matrix  for  a 
conservative  system  is  a  canonical  transformation  in  time  and, 
therefore,  satisfies  the  symplectic  condition.  The  use  of  this 
property  can  then  replace  the  matrix  inversion  in  Eq.  (a.3)  with  a 
simple  rearranging  via  Eq.  (a.13).  If  the  perturbing  force  is  not 


conservative,  as  in  the  case  of  atmospheric  drag,  then  the  <5>(t,t0) 
matrix  must  be  ir verted  numerically. 


A’2  Conversj[on  of  *ttui,t0)  to  »(tt+1,t1) 


To  compute  the  sequential  form  of 

relationships  for  constant  z  , 

ij>  ,  note  the  following 

x^l)  -  4>(t1,t0)x(t0)  ♦  ti)(t1,t0)z0 

(A. 14) 

x(t2)  „  *(t2,t0)x(to)  ♦  n»(t2,t0)z0 

(A. 15) 

x(t2)  -  <J>(t2,t1)x(t1 )  +  ti)(t2,t1)z0 

(A. 16) 

and  substitute  Eq.  (a.14)  fo.'  x ( t ^ )  in  Eq.  (a.Hi) 

x^2)  “  <!*( t -j  )'}>(t -j , tg )x ( t^ )  +  ^ 1 2 •  ) 'i' ( 1 1  •1-0)zo 

+  *(t2.t1)z0 

“  4,(t2.to)5c(to)  *  »t0^zo  +  't’(t2*ti)zo  U*1?) 

Subtract  Eq.  (a.15) 

0  “  )l»(t1ft0)zo  +  'J'(t2,t1)z0  -  ip(t2,t0)z0  (A. 18) 

then  the  equation  for  V'(t2*tj)  is 

-  *(t2.t0)  '  ♦(t2,t1)*(t1,t0)  (A. 19) 


I 


or,  in  general, 


JWU-  JJUIJUUJO 
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■  ^(tj+^,t0)  ~  ^(tj+^  ,tj  )ij/(tj  ,tQ)  (a.20) 

Eqs.  (A. 3)  and  (a.20)  are  used  in  conjunction  with  a 
UTOPIA  run  to  produce  a  file  containing  the  user  spacecraft  state, 
P  and  ijj  matrices  for  each  time  step  desired  during  the  period  of 


interest. 


APPENDIX  B 


THE  U-D  FILTER 

The  GPS  onboard  navigation  filter  is  required  to  perforin 
high-precision  calculations  even  though  it  is  operating  on  a 
microprocessor.  Since  typical  spaceborne  processors  use  small  word 
lengths  (l6  or  32  bits),  techniques  must  be  used  to  improve  the 
numerical  precision  and  stability  of  the  filter.  Several  techniques 
have  been  devised,  all  using  some  form  of  the  square  root  of  the 
covariance.  Maybeck  [ 1 979 ]  and  Tapley,  et  al.  [ 1 980 )  provide 
descriptions  of  the  various  algorithms,  and  the  analysis  by  Tapley, 
et  al.,  shows  that  the  U-D  algorithm  by  Bicrman  [ 1 977 ]  had  the 
lowest  total  numerical  operations,  while  potentially  maintaining 
stability  in  short  word  length  machines.  The  U-D  algorithm  Is 
recommended  for  GPS  navigation  applications,  so  it  was  tested  here 
in  the  autonomous  navigation  role.  The  filter  described  in  this 
report  Includes  model  noise  compensation  (q  matrix)  for  the  reasons 
discussed  in  Section  ^.2.  While  consider  covariance  techniques  are 
appropriate  for  pre-launch  studies,  the  added  complexity  and  core 
storage  required  by  the  consider  parameters  make  it  inefficient  for 
real-time  application. 

As  described  by  Tapley  and  Peters  [ 1 980 ] ,  the  covariance 
matrix  P(tk)  can  be  factored  into 
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Pk.UkDkuJ  (B.l) 

where  U  is  an  upper  triangular  unitary  matrix  (with  ones  on  the 
diagonal)  and  D  is  a  diagonal  matrix.  It  is  well  known  that  this 
factorization  exists  and  is  unique,  even  though  the  matrices  U  and 
D  are  not  unique  [Maybeck,  1979,  p.  392]. 


B.l  U-p  Propagation  via  the  State  Transition  Matrix 

The  matrix  P^  iS  then  propagated  to  time  tk+i  by  the 
state  transition  method  described  in  Section  4.3.  If  the  state 
noise  matrix,  q(t)  ,  is  diagonal  and  is  expressed  in  the  state 
vector  space,  then  b(t)  »  I  and 


k+1 


'k+1 


0(tk+1,tk)Pk^T(tk+1.tk)  +  ♦(tk+1,T)Q(x)^T(tk+1.T)dT 


“  ^tk»1.tk)ukDUy(tk+1,tk) 

Vi 

<f>(tk  +  1 ,  T  )q(  t  )<i>T(  tk+1  ,t)dT  (b.2) 


where 


♦U.tk)  -  A(t)*(t,tk)  ,  ♦(tk,tk)  -  I 


(b.  3 ) 


and 


9FCx(t),t) 


(B.4) 


A  ( t )  = 
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To  propagate  U  and  D  ,  form  the  augmented  matrices 


Wk+1  “  [ lk  +  l  (tk^k  I  ck+l  J 


(b.  5 ) 


and 


Dk  « 


o  Qk/At 


( B .  6 ) 


where 


'k+1 

'k+1  '  J  ^k+l  •T)d't 
t„ 


At  “  fck+1  ‘  bk 


( B  -  7 ) 

( B.  8 ) 


Now  Wk+1  will  not  be  upper  triangular,  but  Uk+1  and  Dk+1  can 
be  obtained  by  a  modified  Gram-Schmidt  orthogonal ization  on  Wk+1 
weighted  by  Kk  such  that 


Pk+1  “  uk+1Dk+1uk+1  “  \  +  1Dkwk+1 


(B.9) 


Note  that  this  method  requires  the  integration  of  the 
n  x  n  i  equation  and  the  n  x  n  Ck+1  integral,  plus  the 
triangularization  of  ¥k+1  . 

B.2  U-D  Propagation  via  U  and  D 

To  improve  the  U-D  propagation  efficiency,  Tapley  and 
Peters  propose  integrating  the  U-D  form  of  the  matrix  Riccati 


4* 


equation 


pU)  -  A(t)F(t)  +  F(t)AT(t)  +  Q (t ) 
If 

P  -  U  F  TJT 
and 

Q  -  Q/2 

then 

-  -1__T  _J._  __±T 

P  ■!)  D  U  +  UDll'I'+UDU 

and 

•  . 

UDIJ1-  t  U  0  ?  «  Atr  +  fr^T  +  jj  + 
Rearranging, 


(U  D  +  U  |  -  Q  U  T  -  A  0  D  )  UT 

*T  1 

+  U  [D  U  +  ^  UT  -  U~ 1  Q  1  -  F  TJT  A 


187 


then  E(t)  +  E^(t)  »  0  and  Eq.  (b. 15)  can  be  simplified  to 

(U  D  +  u|  -  A  ITT?)  UT  -  $  *  E(t)  =  E(t)  (B.  17) 

The  elements  of  E ( t )  can  be  specified  to  maintain  the  triangular 
form  of  U  during  the  integration  by  defining  the  matrices: 

T  2  A  U  D  (B. 18) 

* 

M  =  U  D  +  U  |  -  T  (B. 19) 

Then 

M  T  “  E  •  Q  +  E  (b.20) 

For  U  and  U  to  be  upper  triangular,  there  are 

n(n-l)/2  unknowns  in  the  skew  symmetric  matrix  E.  The  products 
•  « 

U  D  and  U  D  are  upper  triangular  creating  n(n+l)/2  unknowns. 
The  whole  system  (Eq.  ( B . 20 )  then  has  ).  +  n(r^f^  )  =  n  x  n 

unknowns,  the  same  as  P  ,  which  can  be  uniquely  determined.  The 
elements  of  Eq.  (b.20)  are 
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"21  •  * 


el,n-l  *  "en,n-l 

eln  e2n  *  *  *  V 


and  the  solution  to  this  equation  for  M  is  obtained  noting  that 


nn  *  3nn 


(B. 21 ) 


'n-1,n  +  mnn  un-1,n  “  en,n-1 


(B. 22) 


nn,n-1  “  _en,n-1 


(B.23) 


and  that  Equations  (b.22)  and  (b. 23 )  can  be  solved  for  the  two 
unknowns  mn  n-1  and  en(n_i  .  This  process  then  proceeds 
backwards  up  through  the  M  matrix  until  M  and  E  have  been 
determined: 

•  • 

U  and  D  are  then  obtained  from 


"  D  +  U  |  »  H  +  T 


(B. 24 ) 


n  L  -  n  uikdkj 
‘ij  *  -  J,  ulk  <lkj  *  2 


r 

¥• 

t 


i 
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'  “u  *jj  *  •  (B-25) 

(j  **  1 , . . . ,n)(i  «  1  ,  •  • .  ,  j  ) 

» 

For  i  -  j,  Uj^  «  1  ,  so  Uy  -  0  and  Eq.  (B.2ij)  becomes 

^ii  “  2^mii  +  tji^  (*  "  1>--*'n)  (B. 26) 

which  provides  the  solution  for  the  diagonal  matrix  T) 

Since  IF  is  upper  triangular,  terms  for  which  i  >  j  are 
zero,  thus  the  only  case  left  is  for  i  <  j  ,  which  gives 

/  djj  (B.27) 

(i  ■  l,...,j_l)  (j  ■  2, . . . ,n) 

Eqs.  (b.26)  and  (b.27)  are  then  used  in  the  derivative  subroutine  of 

a  numerical  integrator  to  provide  U( t )  and  D(t)  .  Note  that 

T 

UDU  is  never  formed  in  the  propagation  algorithm. 


ij 


mij  +  tij  ‘  uij 


11 

2 


B. 3  U-D  Measurement  Update 

The  Kalman  measurement  update  to  the  covariance  matrix  for 
a  scalar  observation  at  t  -  tk+1  is 


-  P  -  KHP 


(B.28) 


fi  < 
1 


P  -  (i  -  kh)p 
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where 

K  -  F  ht(h  P  HT  ♦  R)-1 
and 

H  -  ^^k+i  ^k+i  ) 

x«(tk+1) 

P  can  be  factored  [Maybeck,  1979]  as 

UDUt  -  TT  F  TJT  -  (l/o)(TT  F  TJT  HT)  H  XT  F  FT 
-  F  [F  -  (l/a)(F  Ft  ht)(dUt  HT)t]  uT  . 

where 

a  -  HPHT  +  R 
Now,  define  the  vectors 

f  -  FT  Ht 

v  -  F  f  -  T5  IJT  Ht 
and  Eq.  (ES.31)  becomes 

UDUt  «  U  [D  -  ( 1  / ct )  v  vT]T  UT 


(B.29) 

(b. 30) 

(B.31) 

(B.32] 

(B. 33) 
(B.3M 

(b. 35) 


where  v  v  and  F  -  (l/o)  v  vT  are  now  symmetric  n  x  n  matrices. 
Applying  the  decomposition  algorithm  to  D  -  (l/o)  v  vT  yields  a 


.MifX 


u  J  1  p.  -  «S® 


...P-PIM-Rid.l1.,' 
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unit  upper  triangular  matrix,  U+  ,  and  a  diagonal  matrix,  D+  . 
Then 

UDUT  «  (U  U+)d  +  (TT  U*)*"  (B.36) 

and 

U  *■  TT  U+  (b. 37 ) 

D  -  D+  (B. 38) 

The  factor  1/a  still  contains  the  product  HPHT  j  however,  by 
calculating  U  and  D  in  a  recursive  manner,  the  explicit 
formation  of  P  is  avoided.  The  method  also  gives  the  Kalman  gain, 
K  ,  in  a  recursive  form.  The  Kalman  state  update  algorithm  then 
proceeds  as  follows: 

Set  a0  “  R  (B. 39) 

Compute 


n 


fi  ■  hi  +  I  hzuii 
£»i+1  x 

(B. 40 ) 

Vi  "  <ii  fi 

(B.41  ) 

“i  ‘  “i-1  +  fi  vi 

(B.«2) 

di  *  *1  °i-1/cti 

(B.U3) 

bi  *  vi 

(b. 44 ) 
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Pi  “  “fi/ai-1 

Uj  ■  u j ^  +  bj  p^  ( J “1 1 2 1 • • • 1 1 - 1  )  (b.^6) 

bj  "  bj  +  UjL  Vi 

The  new  vector  B,  is  then  used  to  compute  the  Kalman  gain 

K  -  B/a  (B.H8) 

n 

and 

8(tk+i)  -x(tk+1)  +  K  [yk+1  -  Hk+1*(tk+1)]  .  (B.49) 

Note  that  this  formulation  still  computes  the  components  of  the 
matrix  P  «  UDU  when  it  is  recursively  solving  for  fj  ,  Vj  and 
.  This  is  true  of  all  square  root  filters,  even  though  the 
product  is  usually  hidden  in  the  recursion  relations. 
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